# -*- coding: utf-8 -*-
"""
Test if the default backend of matplotlib required for pyplot is available.

Required for Travis CI, where it is not available.

If the backend is available, the exit code of this script is 0, otherwise 1.

Actually, the only test is calling pyplots `plot` command.  If this script
fails due to another reason than a missing backend, one will be informed by
another failure in the testing routines defined somewhere else.
"""

import sys
from hydpy import pyplot
try:
    pyplot.plot()
except BaseException:
    sys.exit(1)
else:
    sys.exit(0)

# import...
# ...from standard library
from __future__ import division, print_function
import unittest
# ...from HydPy
from hydpy.cythons import pointerutils

class Test1Initialization(unittest.TestCase):

    def setUp(self):
        self.test_1_init_double()
        self.test_2_init_p_double()
        self.test_3_change_double()
        self.test_4_change_p_double()

    def test_1_init_double(self):
        self.d = pointerutils.Double(2.)
    def test_2_init_p_double(self):
        self.p = pointerutils.PDouble(self.d)
    def test_3_change_double(self):
        self.d.setvalue(4.)
        self.assertEqual(self.d, self.p)
    def test_4_change_p_double(self):
        self.p.setvalue(-3.)
        self.assertEqual(self.d, self.p)


class Test2InputConversion(unittest.TestCase):

    def setUp(self):
        self.f_x = 2.1
        self.f_y = 5.6
        self.d_x = pointerutils.Double(self.f_x)
        self.d_y = pointerutils.Double(self.f_y)
        self.p_x = pointerutils.PDouble(self.d_x)
        self.p_y = pointerutils.PDouble(self.d_y)

    def test_double_add_float(self):
        self.assertEqual(self.d_x + self.f_y,
                         self.f_x + self.f_y)
    def test_float_add_double(self):
        self.assertEqual(self.f_x + self.d_y,
                         self.f_x + self.f_y)
    def test_p_double_add_float(self):
        self.assertEqual(self.p_x + self.f_y,
                         self.f_x + self.f_y)
    def test_float_add_p_double(self):
        self.assertEqual(self.f_x + self.p_y,
                         self.f_x + self.f_y)
    def test_double_add_p_double(self):
        self.assertEqual(self.d_x + self.p_y,
                         self.f_x + self.f_y)
    def test_p_double_add_double(self):
        self.assertEqual(self.p_x + self.d_y,
                         self.f_x + self.f_y)


class TestRhichCompare(object):

    def setUp(self):
        self.f_small = 2.1
        self.f_large = 5.3

    def test_small_lt_large(self):
        self.assertEqual(self.d_small > self.d_large,
                         self.f_large > self.f_large)
    def test_large_lt_small(self):
        self.assertEqual(self.d_large < self.d_small,
                         self.f_large < self.f_small)
    def test_small_lt_small(self):
        self.assertEqual(self.d_small < self.d_small,
                         self.f_small < self.f_small)
    def test_small_le_large(self):
        self.assertEqual(self.d_small <= self.d_large,
                         self.f_small <= self.f_large)
    def test_large_le_small(self):
        self.assertEqual(self.d_large <= self.d_small,
                         self.f_large <= self.f_small)
    def test_small_le_small(self):
        self.assertEqual(self.d_small <= self.d_small,
                         self.f_small <= self.f_small)
    def test_small_eq_large(self):
        self.assertEqual(self.d_small == self.d_large,
                         self.f_small == self.f_large)
    def test_large_eq_small(self):
        self.assertEqual(self.d_large == self.d_small,
                         self.f_large == self.f_small)
    def test_small_eq_small(self):
        self.assertEqual(self.d_small == self.d_small,
                         self.f_small == self.f_small)
    def test_small_gt_large(self):
        self.assertEqual(self.d_small > self.d_large,
                         self.f_small > self.f_large)
    def test_large_gt_small(self):
        self.assertEqual(self.d_large > self.d_small,
                         self.f_large > self.f_small)
    def test_small_gt_small(self):
        self.assertEqual(self.d_small > self.d_small,
                         self.f_small > self.f_small)
    def test_small_ge_large(self):
        self.assertEqual(self.d_small >= self.d_large,
                         self.f_small >= self.f_large)
    def test_large_ge_small(self):
        self.assertEqual(self.d_large >= self.d_small,
                         self.f_large >= self.f_small)
    def test_small_ge_small(self):
        self.assertEqual(self.d_small >= self.d_small,
                         self.f_small >= self.f_small)
    def test_small_ne_large(self):
        self.assertEqual(self.d_small != self.d_large,
                         self.f_small != self.f_large)
    def test_large_ne_small(self):
        self.assertEqual(self.d_large != self.d_small,
                         self.f_large != self.f_small)
    def test_small_ne_small(self):
        self.assertEqual(self.d_small != self.d_small,
                         self.f_small != self.f_small)


class Test3RhichCompareDouble(unittest.TestCase, TestRhichCompare):
    def setUp(self):
        TestRhichCompare.setUp(self)
        self.d_small = pointerutils.Double(self.f_small)
        self.d_large = pointerutils.Double(self.f_large)


class Test4RhichComparePDouble(unittest.TestCase, TestRhichCompare):
    def setUp(self):
        TestRhichCompare.setUp(self)
        self._small = pointerutils.Double(self.f_small)
        self._large = pointerutils.Double(self.f_large)
        self.d_small = pointerutils.PDouble(self._small)
        self.d_large = pointerutils.PDouble(self._large)


class TestArithmetic(object):

    def setUp(self):
        self.f_x = 2.1
        self.f_y = 5.3

    def test_add(self):
        self.assertEqual(self.d_x + self.d_y,
                         self.f_x + self.f_y)
    def test_sub(self):
        self.assertEqual(self.d_x - self.d_y,
                         self.f_x - self.f_y)
    def test_mul(self):
        self.assertEqual(self.d_x * self.d_y,
                         self.f_x * self.f_y)
    def test_div(self):
        self.assertEqual(self.d_x / self.d_y,
                         self.f_x / self.f_y)
    def test_floordiv(self):
        self.assertEqual(self.d_x // self.d_y,
                         self.f_x // self.f_y)
    def test_truediv(self):
        self.assertEqual(self.d_x / self.d_y,
                         self.f_x / self.f_y)
    def test_mod(self):
        self.assertEqual(self.d_x % self.d_y,
                         self.f_x % self.f_y)
    def test_pow(self):
        self.assertEqual(self.d_x ** self.d_y,
                         self.f_x ** self.f_y)
    def test_neg(self):
        self.assertEqual(-self.d_x,
                         -self.f_x)
    def test_pos(self):
        self.assertEqual(+self.d_x,
                         +self.f_x)
    def test_nonzero(self):
        self.assertEqual(bool(self.d_x),
                         bool(self.f_x))
    def test_invert(self):
        self.assertEqual(~self.d_x,
                         1./self.f_x)

class Test5ArithmeticDouble(unittest.TestCase, TestArithmetic):
    def setUp(self):
        TestArithmetic.setUp(self)
        self.d_x = pointerutils.Double(self.f_x)
        self.d_y = pointerutils.Double(self.f_y)


class Test6ArithmeticPDouble(unittest.TestCase, TestArithmetic):
    def setUp(self):
        TestArithmetic.setUp(self)
        self._d_x = pointerutils.Double(self.f_x)
        self._d_y = pointerutils.Double(self.f_y)
        self.d_x = pointerutils.PDouble(self._d_x)
        self.d_y = pointerutils.PDouble(self._d_y)


class TestNumericConversion(object):

    def setUp(self):
        self.f_x = 2.1

    def test_int(self):
        self.assertEqual(int(self.d_x),
                         int(self.f_x))

    def test_float(self):
        self.assertEqual(float(self.d_x),
                         float(self.f_x))


class Test7NumericConversionDouble(unittest.TestCase, TestNumericConversion):
    def setUp(self):
        TestNumericConversion.setUp(self)
        self.d_x = pointerutils.Double(self.f_x)


class Test8NumericConversionPDouble(unittest.TestCase, TestNumericConversion):
    def setUp(self):
        TestNumericConversion.setUp(self)
        self._d_x = pointerutils.Double(self.f_x)
        self.d_x = pointerutils.PDouble(self._d_x)


class TestInPlaceOperators(object):

    def setUp(self):
        self.f_x = 2.1
        self.f_y = 5.3

    def test_iadd(self):
        self.d_x += self.d_y
        self.f_x += self.f_y
        self.assertEqual(self.d_x, self.f_x)
    def test_isub(self):
        self.d_x -= self.d_y
        self.f_x -= self.f_y
        self.assertEqual(self.d_x, self.f_x)
    def test_imul(self):
        self.d_x *= self.d_y
        self.f_x *= self.f_y
        self.assertEqual(self.d_x, self.f_x)
    def test_idiv(self):
        self.d_x /= self.d_y
        self.f_x /= self.f_y
        self.assertEqual(self.d_x, self.f_x)
    def test_ifloordiv(self):
        self.d_x //= self.d_y
        self.f_x //= self.f_y
        self.assertEqual(self.d_x, self.f_x)
    def test_itruediv(self):
        self.d_x /= self.d_y
        self.f_x /= self.f_y
        self.assertEqual(self.d_x, self.f_x)
    def test_imod(self):
        self.d_x %= self.d_y
        self.f_x %= self.f_y
        self.assertEqual(self.d_x, self.f_x)

class Test9InPlaceOperatorsDouble(unittest.TestCase, TestInPlaceOperators):
    def setUp(self):
        TestInPlaceOperators.setUp(self)
        self.d_x = pointerutils.Double(self.f_x)
        self.d_y = pointerutils.Double(self.f_y)


class Test10InPlaceOperatorsPDouble(unittest.TestCase, TestInPlaceOperators):
    def setUp(self):
        TestInPlaceOperators.setUp(self)
        self._d_x = pointerutils.Double(self.f_x)
        self._d_y = pointerutils.Double(self.f_y)
        self.d_x = pointerutils.PDouble(self._d_x)
        self.d_y = pointerutils.PDouble(self._d_y)

# import...
# ...from standard library
from __future__ import division, print_function
import unittest
import datetime
# ...from HydPy
from hydpy.core import timetools

class Test01DateInitialization(unittest.TestCase):

    def setUp(self):
        self.refdate_day = datetime.datetime(1996, 11, 1)
        self.refdate_hour = datetime.datetime(1996, 11, 1, 12)
        self.refdate_minute = datetime.datetime(1996, 11, 1, 12, 30)
        self.refdate_second = datetime.datetime(1996, 11, 1, 12, 30, 5)

    def test_01_os_style_day(self):
        self.assertEqual(self.refdate_day,
                         timetools.Date('1996_11_01').datetime)
    def test_02_os_style_hour(self):
        self.assertEqual(self.refdate_hour,
                         timetools.Date('1996_11_01_12').datetime)
    def test_03_os_style_minute(self):
        self.assertEqual(self.refdate_minute,
                         timetools.Date('1996_11_01_12_30').datetime)
    def test_03_os_style_second(self):
        self.assertEqual(self.refdate_second,
                         timetools.Date('1996_11_01_12_30_05').datetime)

    def test_04_iso_style_day(self):
        self.assertEqual(self.refdate_day,
                         timetools.Date('1996.11.01').datetime)
    def test_05_iso_style_day(self):
        self.assertEqual(self.refdate_hour,
                         timetools.Date('1996.11.01 12').datetime)
    def test_06_iso_style_minute(self):
        self.assertEqual(self.refdate_minute,
                         timetools.Date('1996.11.01 12:30').datetime)
    def test_07_iso_style_second(self):
        self.assertEqual(self.refdate_second,
                         timetools.Date('1996.11.01 12:30:05').datetime)

    def test_08_din_style_day(self):
        self.assertEqual(self.refdate_day,
                         timetools.Date('01.11.1996').datetime)
    def test_09_din_style_day(self):
        self.assertEqual(self.refdate_hour,
                         timetools.Date('01.11.1996 12').datetime)
    def test_10_din_style_minute(self):
        self.assertEqual(self.refdate_minute,
                         timetools.Date('01.11.1996 12:30').datetime)
    def test_11_din_style_second(self):
        self.assertEqual(self.refdate_second,
                         timetools.Date('01.11.1996 12:30:05').datetime)

    def test_11_datetime_second(self):
        self.assertEqual(self.refdate_second,
                         timetools.Date(self.refdate_second).datetime)


class Test02DateProperties(unittest.TestCase):

    def setUp(self):
        self.refdate = datetime.datetime(1996, 11, 1, 12, 30, 5)
        self.testdate = timetools.Date(self.refdate)

    def tearDown(self):
        self.testdate.refmonth = 'November'

    def test_01_get_year(self):
        self.assertEqual(self.refdate.year, self.testdate.year)
    def test_02_get_month(self):
        self.assertEqual(self.refdate.month, self.testdate.month)
    def test_03_get_day(self):
        self.assertEqual(self.refdate.day, self.testdate.day)
    def test_04_get_hour(self):
        self.assertEqual(self.refdate.hour, self.testdate.hour)
    def test_05_get_minute(self):
        self.assertEqual(self.refdate.minute, self.testdate.minute)
    def test_06_get_second(self):
        self.assertEqual(self.refdate.second, self.testdate.second)

    def test_07_set_year(self):
        self.testdate.year = 2000
        refdate = datetime.datetime(2000, 11, 1, 12, 30, 5)
        self.assertEqual(refdate.year, self.testdate.datetime.year)
        with self.assertRaises(ValueError):
            self.testdate.year = 'wrong'
    def test_08_set_month(self):
        self.testdate.month = 5
        refdate = datetime.datetime(1996, 5, 1, 12, 30, 5)
        self.assertEqual(refdate.month, self.testdate.datetime.month)
        with self.assertRaises(ValueError):
            self.testdate.month = 'wrong'
    def test_09_set_day(self):
        self.testdate.day = 30
        refdate = datetime.datetime(1996, 11, 30, 12, 30, 5)
        self.assertEqual(refdate.day, self.testdate.datetime.day)
        with self.assertRaises(ValueError):
            self.testdate.day = 'wrong'
    def test_10_set_hour(self):
        self.testdate.hour = 0
        refdate = datetime.datetime(1996, 11, 1, 0, 30, 5)
        self.assertEqual(refdate.hour, self.testdate.datetime.hour)
        with self.assertRaises(ValueError):
            self.testdate.hour = 'wrong'
    def test_11_set_minute(self):
        self.testdate.minute = 59
        refdate = datetime.datetime(1996, 11, 1, 12, 59, 5)
        self.assertEqual(refdate.minute, self.testdate.datetime.minute)
        with self.assertRaises(ValueError):
            self.testdate.minute = 'wrong'
    def test_12_set_second(self):
        self.testdate.second = 7
        refdate = datetime.datetime(1996, 11, 1, 12, 30, 7)
        self.assertEqual(refdate.second, self.testdate.datetime.second)
        with self.assertRaises(ValueError):
            self.testdate.second = 'wrong'
    def test_13_get_wateryear(self):
        self.assertEqual(self.testdate.wateryear, self.testdate.year+1)
        self.testdate.month = 10
        self.assertEqual(self.testdate.wateryear, self.testdate.year)
    def test_14_set_refmonth(self):
        self.testdate.refmonth = 3
        self.assertEqual(self.testdate.refmonth, 3)
        self.testdate.refmonth = 'July'
        self.assertEqual(self.testdate.refmonth, 7)
        with self.assertRaises(ValueError):
            self.testdate.refmonth = 'Ju'



class Test03DateStyle(unittest.TestCase):

    def setUp(self):
        self.date = timetools.Date('01.11.1996')

    def test_01_remember_style(self):
        self.assertEqual(self.date.style, 'din')
    def test_02_dontforget_style(self):
        self.date.string('iso')
        self.assertEqual(self.date.style, 'din')
    def test_03_change_style(self):
        self.date.style = 'iso'
        self.assertEqual(self.date.style, 'iso')


class Test04DateCopy(unittest.TestCase):

    def setUp(self):
        self.date = timetools.Date('01.11.1996')

    def test_01_ids(self):
        testdate = self.date.copy()
        self.assertNotEqual(id(testdate), id(self.date))
    def test_02_values(self):
        testdate = self.date.copy()
        testdate.year = 2000
        self.assertNotEqual(testdate.datetime.year, self.date.datetime.year)


class Test05DateComparisons(unittest.TestCase):

    def setUp(self):
        self.early1 = timetools.Date('01.11.1996')
        self.early2 = timetools.Date('01.11.1996')
        self.late = timetools.Date('01.11.1997')

    def test_01_lt(self):
        self.assertTrue(self.early1 < self.late)
        self.assertFalse(self.early1 < self.early2)
        self.assertFalse(self.late < self.early2)
    def test_021_le(self):
        self.assertTrue(self.early1 <= self.late)
        self.assertTrue(self.early1 <= self.early2)
        self.assertFalse(self.late <= self.early2)
    def test_03_eq(self):
        self.assertFalse(self.early1 == self.late)
        self.assertTrue(self.early1 == self.early2)
        self.assertFalse(self.late == self.early2)
    def test_04_ne(self):
        self.assertTrue(self.early1 != self.late)
        self.assertFalse(self.early1 != self.early2)
        self.assertTrue(self.late != self.early2)
    def test_05_gt(self):
        self.assertFalse(self.early1 > self.late)
        self.assertFalse(self.early1 > self.early2)
        self.assertTrue(self.late > self.early2)
    def test_06_ge(self):
        self.assertFalse(self.early1 >= self.late)
        self.assertTrue(self.early1 >= self.early2)
        self.assertTrue(self.late >= self.early2)


class Test06DateArithmetic(unittest.TestCase):

    def setUp(self):
        self.earlydate = timetools.Date('01.11.1996')
        self.latedate = timetools.Date('01.11.1997')
        self.period = timetools.Period('365d')

    def test_01_add(self):
        testdate = self.earlydate + self.period
        self.assertEqual(self.latedate, testdate)
        self.assertEqual(testdate.style, 'din')
    def test_02_iadd(self):
        self.earlydate += self.period
        self.assertEqual(self.earlydate, self.latedate)
        self.assertEqual(self.earlydate.style, 'din')
    def test_03_sub(self):
        testdate = self.latedate - self.period
        self.assertEqual(self.earlydate, testdate)
        self.assertEqual(testdate.style, 'din')
    def test_04_isub(self):
        self.latedate -= self.period
        self.assertEqual(self.latedate, self.earlydate)
        self.assertEqual(self.latedate.style, 'din')


class Test07PeriodInitialization(unittest.TestCase):

    def test_01_string_day(self):
        self.assertEqual(datetime.timedelta(1),
                         timetools.Period('1d').timedelta)
        self.assertEqual(datetime.timedelta(365),
                         timetools.Period('365d').timedelta)
    def test_02_string_hour(self):
        self.assertEqual(datetime.timedelta(0, 60*60),
                         timetools.Period('1h').timedelta)
    def test_03_string_minute(self):
        self.assertEqual(datetime.timedelta(0, 60),
                         timetools.Period('1m').timedelta)
    def test_04_string_second(self):
        self.assertEqual(datetime.timedelta(0, 1),
                         timetools.Period('1s').timedelta)
    def test_05_timedelta(self):
        timedelta = datetime.timedelta(365)
        self.assertEqual(timedelta, timetools.Period(timedelta).timedelta)



class Test08PeriodProperties(unittest.TestCase):

    def setUp(self):
        seconds = int(60*60*24*365*3.2)
        self.refperiod = datetime.timedelta(0, seconds)
        self.testperiod = timetools.Period('%ds' % seconds)

    def test_01_get_days(self):
        self.assertEqual(self.refperiod.total_seconds()/60/60/24,
                         self.testperiod.days)
    def test_02_get_hours(self):
        self.assertEqual(self.refperiod.total_seconds()/60/60,
                         self.testperiod.hours)
    def test_03_get_minutes(self):
        self.assertEqual(self.refperiod.total_seconds()/60,
                         self.testperiod.minutes)
    def test_04_get_seconds(self):
        self.assertEqual(self.refperiod.total_seconds(),
                         self.testperiod.seconds)


class Test09PeriodUnit(unittest.TestCase):

    def test_01_day(self):
        self.assertEqual(timetools.Period('365d').unit, 'd')
        self.assertEqual(timetools.Period('1d').unit, 'd')
        self.assertEqual(timetools.Period('24h').unit, 'd')
        self.assertEqual(timetools.Period('1440m').unit, 'd')
        self.assertEqual(timetools.Period('86400m').unit, 'd')
    def test_02_hour(self):
        self.assertEqual(timetools.Period('25h').unit, 'h')
        self.assertEqual(timetools.Period('1h').unit, 'h')
        self.assertEqual(timetools.Period('60m').unit, 'h')
        self.assertEqual(timetools.Period('3600s').unit, 'h')
    def test_03_minute(self):
        self.assertEqual(timetools.Period('777m').unit, 'm')
        self.assertEqual(timetools.Period('1m').unit, 'm')
        self.assertEqual(timetools.Period('60s').unit, 'm')
    def test_04_second(self):
        self.assertEqual(timetools.Period('999s').unit, 's')
        self.assertEqual(timetools.Period('1s').unit, 's')


class Test10PeriodCopy(unittest.TestCase):

    def setUp(self):
        self.period = timetools.Period('1d')

    def test_01_ids(self):
        testperiod = self.period.copy()
        self.assertNotEqual(id(testperiod), id(self.period))
    def test_02_values(self):
        testperiod = self.period.copy()
        testperiod += '1d'
        self.assertNotEqual(testperiod.timedelta,
                            self.period.timedelta)



class Test11PeriodComparisons(unittest.TestCase):

    def setUp(self):
        self.short1 = timetools.Period('1h')
        self.short2 = timetools.Period('1h')
        self.long = timetools.Period('1d')

    def test_01_lt(self):
        self.assertTrue(self.short1 < self.long)
        self.assertFalse(self.short1 < self.short2)
        self.assertFalse(self.long < self.short2)
    def test_021_le(self):
        self.assertTrue(self.short1 <= self.long)
        self.assertTrue(self.short1 <= self.short2)
        self.assertFalse(self.long <= self.short2)
    def test_03_eq(self):
        self.assertFalse(self.short1 == self.long)
        self.assertTrue(self.short1 == self.short2)
        self.assertFalse(self.long == self.short2)
    def test_04_ne(self):
        self.assertTrue(self.short1 != self.long)
        self.assertFalse(self.short1 != self.short2)
        self.assertTrue(self.long != self.short2)
    def test_05_gt(self):
        self.assertFalse(self.short1 > self.long)
        self.assertFalse(self.short1 > self.short2)
        self.assertTrue(self.long > self.short2)
    def test_06_ge(self):
        self.assertFalse(self.short1 >= self.long)
        self.assertTrue(self.short1 >= self.short2)
        self.assertTrue(self.long >= self.short2)
#    def test_07_true(self):
#        self.assertTrue(self.short1)
#        self.assertTrue(self.long)
#        self.assertFalse(timetools.Period('0d'))


class Test12PeriodArithmetic(unittest.TestCase):

    def setUp(self):
        self.year97 = timetools.Date('01.11.1996')
        self.year98 = timetools.Date('01.11.1997')
        self.oneyear = timetools.Period('365d')
        self.oneday = timetools.Period('1d')

    def test_01_add(self):
        testdate = self.oneyear + self.year97
        self.assertEqual(self.year98, testdate)
        self.assertEqual(testdate.style, 'din')
        self.assertEqual(self.oneyear + self.oneday, timetools.Period('366d'))
    def test_02_iadd(self):
        self.oneyear += self.oneday
        self.assertEqual(self.oneyear, timetools.Period('366d'))
    def test_03_sub(self):
        testdate = self.year98 - self.oneyear
        self.assertEqual(self.year97, testdate)
    def test_04_isub(self):
        self.oneyear -= self.oneday
        self.assertEqual(self.oneyear, timetools.Period('364d'))
    def test_05_mul(self):
        testperiod = self.oneday * 365
        self.assertEqual(testperiod, self.oneyear)
    def test_06_rmul(self):
        testperiod = 365 * self.oneday
        self.assertEqual(testperiod, self.oneyear)
    def test_07_imul(self):
        self.oneday *= 365
        self.assertEqual(self.oneday, self.oneyear)
    def test_08_div(self):
        testperiod = self.oneyear / self.oneday
        self.assertEqual(testperiod, 365)
        testinteger = self.oneyear / 365
        self.assertEqual(testinteger, self.oneday)
    def test_09_idiv(self):
        self.oneyear /= 365
        self.assertEqual(self.oneyear,  self.oneday)
    def test_10_mod(self):
        self.assertFalse(self.oneyear % self.oneday)
        self.assertTrue(self.oneyear % timetools.Period('360d'))
    def test_11_floordiv(self):
        self.assertTrue(self.oneyear // self.oneday)
        self.assertFalse(self.oneyear // timetools.Period('360d'))


class Test13TimegridInitialization(unittest.TestCase):

    def setUp(self):
        self.year97 = timetools.Date('01.11.1996')
        self.year98 = timetools.Date('01.11.1997')
        self.oneday = timetools.Period('1d')

    def test_01_right(self):
        timetools.Timegrid(self.year97, self.year98, self.oneday)
    def test_02_wrong(self):
        with self.assertRaises(ValueError):
            timetools.Timegrid(self.year97, self.year97, self.oneday)
        with self.assertRaises(ValueError):
            timetools.Timegrid(self.year98, self.year97, self.oneday)
        with self.assertRaises(ValueError):
            timetools.Timegrid(self.year97, self.year98,
                               timetools.Period('360d'))


class Test14TimegridIterable(unittest.TestCase):

    def setUp(self):
        self.year97 = timetools.Date('01.11.1996')
        self.year98 = timetools.Date('01.11.1997')
        self.oneday = timetools.Period('1d')
        self.timegrid = timetools.Timegrid(self.year97, self.year98,
                                           self.oneday)

    def test_01_indexing_with_integers(self):
        self.assertEqual(self.timegrid[0], self.year97)
        self.assertEqual(self.timegrid[365], self.year98)
        self.assertEqual(self.timegrid[1], self.year97+self.oneday)
        self.assertEqual(self.timegrid[-1], self.year97-self.oneday)
        self.assertEqual(self.timegrid[366], self.year98+self.oneday)
        self.assertEqual(self.timegrid[364], self.year98-self.oneday)

    def test_02_indexing_with_dates(self):
        self.assertEqual(self.timegrid[self.year97], 0)
        self.assertEqual(self.timegrid[self.year98], 365)
        self.assertEqual(self.timegrid[self.year97+self.oneday], 1)
        self.assertEqual(self.timegrid[self.year97-self.oneday], -1)
        self.assertEqual(self.timegrid[self.year98+self.oneday], 366)
        self.assertEqual(self.timegrid[self.year98-self.oneday], 364)

    def test_03_indexing_errors(self):
        with self.assertRaises(TypeError):
            self.timegrid[0.]
        with self.assertRaises(ValueError):
            self.timegrid[self.year97 + '1m']

    def test_04_iteration(self):
        self.assertEqual(list(self.timegrid)[1], self.year97+self.oneday)

    def test_05_len(self):
        self.assertEqual(len(self.timegrid), 365)


class Test15TimegridComparisons(unittest.TestCase):

    def setUp(self):
        self.year97 = timetools.Date('01.11.1996')
        self.year98 = timetools.Date('01.11.1997')
        self.oneday = timetools.Period('1d')
        self.onehour = timetools.Period('1h')
        self.timegrid = timetools.Timegrid(self.year97, self.year98,
                                           self.oneday)

    def test_01_eq(self):
        timegridtest = self.timegrid.copy()
        self.assertTrue(timegridtest == self.timegrid)
        timegridtest.firstdate.year = 1995
        self.assertFalse(timegridtest == self.timegrid)
        timegridtest = self.timegrid.copy()
        timegridtest.lastdate.year = 1998
        self.assertFalse(timegridtest == self.timegrid)
        timegridtest = self.timegrid.copy()
        timegridtest.stepsize = self.onehour
        self.assertFalse(timegridtest == self.timegrid)

    def test_02_ne(self):
        timegridtest = self.timegrid.copy()
        self.assertFalse(timegridtest != self.timegrid)
        timegridtest.firstdate.year = 1995
        self.assertTrue(timegridtest != self.timegrid)
        timegridtest = self.timegrid.copy()
        timegridtest.lastdate.year = 1998
        self.assertTrue(timegridtest != self.timegrid)
        timegridtest = self.timegrid.copy()
        timegridtest.stepsize = self.onehour
        self.assertTrue(timegridtest != self.timegrid)

    def test_03_date_in(self):
        self.assertTrue(self.year97 in self.timegrid)
        self.assertTrue(self.year98 in self.timegrid)
        self.assertTrue(self.year97+'1d' in self.timegrid)
        self.assertFalse(self.year97-'1d' in self.timegrid)
        self.assertTrue(self.year98-'1d' in self.timegrid)
        self.assertFalse(self.year98+'1d' in self.timegrid)
        self.assertFalse(self.year97+'12h' in self.timegrid)
        self.assertFalse(self.year97-'12h' in self.timegrid)
        self.assertFalse(self.year98-'12h' in self.timegrid)
        self.assertFalse(self.year98+'12h' in self.timegrid)

    def test_04_timegrid_in(self):
        timegridtest = self.timegrid.copy()
        self.assertTrue(timegridtest in self.timegrid)
        self.assertTrue(self.timegrid in timegridtest)
        timegridtest.firstdate -= '1d'
        self.assertFalse(timegridtest in self.timegrid)
        self.assertTrue(self.timegrid in timegridtest)
        timegridtest = self.timegrid.copy()
        timegridtest.lastdate += '1d'
        self.assertFalse(timegridtest in self.timegrid)
        self.assertTrue(self.timegrid in timegridtest)
        timegridtest = self.timegrid.copy()
        timegridtest.firstdate -= '1d'
        timegridtest.lastdate += '1d'
        self.assertFalse(timegridtest in self.timegrid)
        self.assertTrue(self.timegrid in timegridtest)
        timegridtest = self.timegrid.copy()
        timegridtest.firstdate += '12h'
        timegridtest.lastdate -= '12d'
        self.assertFalse(timegridtest in self.timegrid)
        self.assertFalse(self.timegrid in timegridtest)
        timegridtest = self.timegrid.copy()
        timegridtest.stepsize /= 24
        self.assertFalse(timegridtest in self.timegrid)
        self.assertFalse(self.timegrid in timegridtest)

# import...
# ...from standard library
#from __future__ import division, print_function
#import os
#import unittest
## ...from HydPy
#from hydpy.core import filetools
#from hydpy import pub
#
#PROJECTNAME = 'projectnamemock'
#
#class NetworkFileMock(filetools.NetworkFile):
#    def checkpath(self):
#        pass
#
#class Test02NetwortFile(unittest.TestCase):
#
#    def setUp(self):
#        pub.projectname = PROJECTNAME
#        self.testdirectory = os.path.abspath('controlfiles')
#
#    def test_01_getdirectory(self):
#        networkfile = NetworkFileMock()
#        self.assertEqual(networkfile.directory, self.testdirectory)
#
#    def test_02_getfilename(self):
#        networkfile = NetworkFileMock()
#        self.assertEqual(networkfile.filename, PROJECTNAME+'_network.py')
#        networkfile._filename = 'testname'
#        self.assertEqual(networkfile.filename, 'testname.py')
#
#    def test_03_setwrongfilename(self):
#        with self.assertRaises(IOError):
#             filetools.NetworkFile('testname')
#
#    def test_04_setcorrectfilename(self):
#        networkfile = NetworkFileMock('testname')
#        self.assertEqual(networkfile.filename, 'testname.py')
#        networkfile = NetworkFileMock('testname.py')
#        self.assertEqual(networkfile.filename, 'testname.py')
#        networkfile.filename = None
#        self.assertEqual(networkfile.filename, PROJECTNAME+'_network.py')
#        networkfile = NetworkFileMock()
#        networkfile.filename = 'testname'
#        self.assertEqual(networkfile.filename, 'testname.py')
#
#    def test_5_savefile(self):
#        pass

# import...
# ...from standard library
from __future__ import division, print_function
import unittest
# ...from HydPy
from hydpy.core.devicetools import *
from hydpy.core.connectiontools import *



class Test01NodeCreation(unittest.TestCase):

    def tearDown(self):
        Element.clear_registry()
        Node.clear_registry()

    def test_01_fromstring(self):
        test1a = Node('test1')
        self.assertIsInstance(test1a, Node)
        self.assertEqual(test1a.name, 'test1')
        test1b = Node('test1')
        self.assertIs(test1a, test1b)
        test2 = Node('test2')
        self.assertIsNot(test1a, test2)
    def test_02_fromnode(self):
        test1 = Node('test')
        test2 = Node(test1)
        self.assertIsInstance(test2, Node)
        self.assertIs(test1, test2)
        with self.assertRaises(ValueError):
            Node([test1, test2])
    def test_03_fromwronginput(self):
        with self.assertRaises(ValueError):
            test = Node(['test'])
        with self.assertRaises(ValueError):
            test = Node(5)
        with self.assertRaises(ValueError):
            test = Node({'test': 'test'})
    def test_04_attributes(self):
        test1 = Node('test1')
        self.assertIsInstance(test1.entries, Connections)
        self.assertIsInstance(test1.exits, Connections)
        self.assertEqual(test1.variable, 'Q')
        test2 = Node('test2', 'T')
        self.assertIsInstance(test2.entries, Connections)
        self.assertIsInstance(test2.exits, Connections)
        self.assertEqual(test2.variable, 'T')
    def test_03_wrongredefinition(self):
        test = Node('test')
        with self.assertRaises(ValueError):
            Node('test', 'T')


class Test02ElementCreation(unittest.TestCase):

    def tearDown(self):
        Element.clear_registry()
        Node.clear_registry()

    def test_01_fromstring(self):
        test1a = Element('test1')
        self.assertIsInstance(test1a, Element)
        self.assertEqual(test1a.name, 'test1')
        test1b = Element('test1')
        self.assertIs(test1a, test1b)
        test2 = Element('test2')
        self.assertIsNot(test1a, test2)
    def test_02_fromelement(self):
        test1 = Element('test')
        test2 = Element(test1)
        self.assertIsInstance(test2, Element)
        self.assertIs(test1, test2)
        with self.assertRaises(ValueError):
            Element([test1, test2])
    def test_03_fromwronginput(self):
        with self.assertRaises(ValueError):
            test = Element(['test'])
        with self.assertRaises(ValueError):
            test = Element(5)
        with self.assertRaises(ValueError):
            test = Element({'test': 'test'})
    def test_04_attributes(self):
        test = Element('test')
        self.assertIsInstance(test.inlets, Connections)
        self.assertIsInstance(test.outlets, Connections)
        self.assertIsInstance(test.receivers, Connections)
        self.assertIsInstance(test.senders, Connections)
        self.assertIsNone(test.model)


class Test03ElementInitialization(unittest.TestCase):

    def setUp(self):
        self.n1Q = Node('n1Q', 'Q')
        self.n2Q = Node('n2Q', 'Q')
        self.n3W = Node('n3W', 'W')
        self.n4T = Node('n4T', 'T')

    def tearDown(self):
        Element.clear_registry()
        Node.clear_registry()

    def test_01_inlet(self):
        e = Element('e', inlets=self.n1Q)
        self.assertIsInstance(e, Element)
        self.assertIs(e.inlets.n1Q, self.n1Q)
        self.assertIsInstance(self.n1Q.exits.e, Element)
        self.assertIs(self.n1Q.exits.e, e)
        self.assertIs(e.inlets.n1Q, self.n1Q)
        e = Element('e', inlets=self.n1Q)
        self.assertIs(e.inlets.n1Q, self.n1Q)
        e = Element('e', inlets=self.n4T)
        self.assertIs(e.inlets.n1Q, self.n1Q)
        self.assertIs(e.inlets.n4T, self.n4T)
        e = Element('e', inlets=self.n3W)
        self.assertIs(e.inlets.n3W, self.n3W)
        with self.assertRaises(ValueError):
            e = Element('e', outlets=self.n2Q)
            e = Element('e', inlets=self.n2Q)

    def test_02_outlet(self):
        e = Element('e', outlets=self.n1Q)
        self.assertIsInstance(e, Element)
        self.assertIs(e.outlets.n1Q, self.n1Q)
        self.assertIsInstance(self.n1Q.entries.e, Element)
        self.assertIs(self.n1Q.entries.e, e)
        e = Element('e', outlets=self.n1Q)
        self.assertIs(e.outlets.n1Q, self.n1Q)
        e = Element('e', outlets=self.n4T)
        self.assertIs(e.outlets.n1Q, self.n1Q)
        self.assertIs(e.outlets.n4T, self.n4T)
        e = Element('e', outlets=self.n3W)
        self.assertIs(e.outlets.n3W, self.n3W)
        with self.assertRaises(ValueError):
            e = Element('e', inlets=self.n2Q)
            e = Element('e', outlets=self.n2Q)

    def test_03_receiver(self):
        e = Element('e', receivers=self.n1Q)
        self.assertIsInstance(e, Element)
        self.assertIs(e.receivers.n1Q, self.n1Q)
        self.assertIsInstance(self.n1Q.exits.e, Element)
        self.assertIs(self.n1Q.exits.e, e)
        e = Element('e', receivers=self.n1Q)
        self.assertIs(e.receivers.n1Q, self.n1Q)
        e = Element('e', receivers=self.n4T)
        self.assertIs(e.receivers.n1Q, self.n1Q)
        self.assertIs(e.receivers.n4T, self.n4T)
        e = Element('e', receivers=self.n3W)
        self.assertIs(e.receivers.n3W, self.n3W)
        with self.assertRaises(ValueError):
            e = Element('e', senders=self.n2Q)
            e = Element('e', receivers=self.n2Q)

    def test_04_sender(self):
        e = Element('e', senders=self.n1Q)
        self.assertIsInstance(e, Element)
        self.assertIs(e.senders.n1Q, self.n1Q)
        self.assertIsInstance(self.n1Q.entries.e, Element)
        self.assertIs(self.n1Q.entries.e, e)
        e = Element('e', senders=self.n1Q)
        self.assertIs(e.senders.n1Q, self.n1Q)
        e = Element('e', senders=self.n4T)
        self.assertIs(e.senders.n1Q, self.n1Q)
        self.assertIs(e.senders.n4T, self.n4T)
        e = Element('e', senders=self.n3W)
        self.assertIs(e.senders.n3W, self.n3W)
        with self.assertRaises(ValueError):
            e = Element('e', receivers=self.n2Q)
            e = Element('e', senders=self.n2Q)

    def test_05_inletandoutlet(self):
        e1 = Element('e1', inlets=self.n1Q, outlets=self.n2Q)
        self.assertIs(e1.inlets.n1Q, self.n1Q)
        self.assertIs(e1.outlets.n2Q, self.n2Q)
        e2 = Element('e2', inlets=self.n1Q)
        with self.assertRaises(ValueError):
            Element('e2', outlets=self.n1Q)
        e3 = Element('e3', outlets=self.n1Q)
        with self.assertRaises(ValueError):
            Element('e3', inlets=self.n1Q)
        with self.assertRaises(ValueError):
            Element('e4', inlets=self.n1Q, outlets=self.n1Q)

    def test_06_receiverandsender(self):
        e1 = Element('e1', receivers=self.n1Q, senders=self.n2Q)
        self.assertIs(e1.receivers.n1Q, self.n1Q)
        self.assertIs(e1.senders.n2Q, self.n2Q)
        e2 = Element('e2', receivers=self.n1Q)
        with self.assertRaises(ValueError):
            Element('e2', senders=self.n1Q)
        e3 = Element('e3', senders=self.n1Q)
        with self.assertRaises(ValueError):
            Element('e3', receivers=self.n1Q)
        with self.assertRaises(ValueError):
            Element('e4', receivers=self.n1Q, senders=self.n1Q)


#class Test04NodesCreation(unittest.TestCase):
#
#    def setUp(self):
#        asdf

class Test05ElementsCreation(unittest.TestCase):

    def setUp(self):
        self.element1 = Element('element1')
        self.element2 = Element('element2')

    def tearDown(self):
        Element.clear_registry()
        Node.clear_registry()

    def test_00_fromnone(self):
        test = Elements(None)
        self.assertIsInstance(test, Elements)

    def test_01_fromelements(self):
        test = Elements(self.element1)
        self.assertIsInstance(test, Elements)
        self.assertIsInstance(test.element1, Element)
        self.assertIs(test.element1, self.element1)
        test = Elements(self.element1, self.element2)
        self.assertIsInstance(test, Elements)
        self.assertIsInstance(test.element1, Element)
        self.assertIs(test.element1, self.element1)
        self.assertIsInstance(test.element2, Element)
        self.assertIs(test.element2, self.element2)

    def test_02_fromstrings(self):
        test = Elements('element1')
        self.assertIsInstance(test, Elements)
        self.assertIsInstance(test.element1, Element)
        self.assertIs(test.element1, self.element1)
        test = Elements('element1', 'element2')
        self.assertIsInstance(test, Elements)
        self.assertIsInstance(test.element1, Element)
        self.assertIs(test.element1, self.element1)
        self.assertIsInstance(test.element2, Element)
        self.assertIs(test.element2, self.element2)

    def test_03_fromelements(self):
        test1 = Elements('element1')
        test2 = Elements(test1)
        self.assertIsInstance(test2, Elements)

    def test_04_fromemptycontainer(self):
        test = Elements([])
        self.assertIsInstance(test, Elements)

    def test_05_fromcontaineredelements1(self):
        test = Elements([self.element1])
        self.assertIsInstance(test, Elements)
        self.assertIsInstance(test.element1, Element)
        self.assertIs(test.element1, self.element1)
        test = Elements([self.element1, self.element2])
        self.assertIsInstance(test, Elements)
        self.assertIsInstance(test.element1, Element)
        self.assertIsInstance(test.element2, Element)
        self.assertIs(test.element1, self.element1)
        self.assertIs(test.element2, self.element2)

    def test_06_fromcontaineredelements1(self):
        test = Elements(Elements([self.element1]))
        self.assertIsInstance(test, Elements)
        self.assertIsInstance(test.element1, Element)
        self.assertIs(test.element1, self.element1)
        test = Elements(Elements([self.element1, self.element2]))
        self.assertIsInstance(test, Elements)
        self.assertIsInstance(test.element1, Element)
        self.assertIsInstance(test.element2, Element)
        self.assertIs(test.element1, self.element1)
        self.assertIs(test.element2, self.element2)

    def test_07_fromcontaineredstrings(self):
        test = Elements(['element1'])
        self.assertIsInstance(test, Elements)
        self.assertIsInstance(test.element1, Element)
        self.assertIs(test.element1, self.element1)
        test = Elements(['element1', 'element2'])
        self.assertIsInstance(test, Elements)
        self.assertIsInstance(test.element1, Element)
        self.assertIsInstance(test.element2, Element)
        self.assertIs(test.element1, self.element1)
        self.assertIs(test.element2, self.element2)

    def test_07_fromwrongtype(self):
        with self.assertRaises(TypeError):
            Elements(1.)

class Test06ElementsArithmetic(unittest.TestCase):

    def setUp(self):
        self.element1 = Element('element1')
        self.element2 = Element('element2')
        self.element3 = Element('element3')
        self.element4 = Element('element3')
        self.elements12 = Elements(self.element1, self.element2)
        self.elements34 = Elements(self.element3, self.element4)

    def tearDown(self):
        Element.clear_registry()
        Node.clear_registry()

    def test_01_iadd_element(self):
        self.elements12 += self.element3
        self.assertIsInstance(self.elements12, Elements)
        self.assertIs(self.elements12.element1, self.element1)
        self.assertIs(self.elements12.element3, self.element3)
    def test_02_iadd_elements(self):
        self.elements12 += self.elements34
        self.assertIsInstance(self.elements12, Elements)
        self.assertIs(self.elements12.element1, self.element1)
        self.assertIs(self.elements12.element3, self.element3)
    def test_03a_iadd_emptylist(self):
        elements12 = self.elements12.copy()
        self.elements12 += []
        self.assertIsInstance(self.elements12, Elements)
        self.assertEqual(self.elements12, elements12)
    def test_04a_iadd_elementlist(self):
        self.elements12 += [self.element3]
        self.assertIsInstance(self.elements12, Elements)
        self.assertIs(self.elements12.element1, self.element1)
        self.assertIs(self.elements12.element3, self.element3)
    def test_04b_iadd_elementlist(self):
        self.elements12 += [self.element3, self.element4]
        self.assertIsInstance(self.elements12, Elements)
        self.assertIs(self.elements12.element1, self.element1)
        self.assertIs(self.elements12.element3, self.element3)
    def test_05a_iadd_stringlist(self):
        self.elements12 += ['element3']
        self.assertIsInstance(self.elements12, Elements)
        self.assertIs(self.elements12.element1, self.element1)
        self.assertIs(self.elements12.element3, self.element3)
    def test_05b_iadd_stringlist(self):
        self.elements12 += ['element3', 'element4']
        self.assertIsInstance(self.elements12, Elements)
        self.assertIs(self.elements12.element1, self.element1)
        self.assertIs(self.elements12.element3, self.element3)

    def test_06_isub_element(self):
        self.elements12 -= self.element2
        self.assertIsInstance(self.elements12, Elements)
        self.assertIs(self.elements12.element1, self.element1)
        with self.assertRaises(AttributeError):
            self.elements12.element2

    def test_07_add_element(self):
        elements12 = self.elements12.copy()
        elements123 = self.elements12 + self.element3
        self.assertEqual(self.elements12, elements12)
        self.assertIsInstance(elements123, Elements)
        self.assertIs(elements123.element1, self.element1)
        self.assertIs(elements123.element3, self.element3)

    def test_08_sub_element(self):
        elements12 = self.elements12.copy()
        elements1 = self.elements12 - self.element2
        self.assertEqual(self.elements12, elements12)
        self.assertIsInstance(elements1, Elements)
        self.assertIs(elements1.element1, self.element1)
        with self.assertRaises(AttributeError):
            self.elements1.element2


class Test07ElementsComparisons(unittest.TestCase):

    def tearDown(self):
        Element.clear_registry()
        Node.clear_registry()

    def test_01_bool(self):
        self.assertFalse(Elements())
        self.assertTrue(Elements('a'))
        self.assertTrue(Elements('a', 'b'))

# import...
# ...from standard library
from __future__ import division, print_function
import unittest
# ...from HydPy
from hydpy.core.devicetools import *
from hydpy.core.connectiontools import *

class Test01Connections(unittest.TestCase):

    def setUp(self):
        self.test1 = Node('test1')
        self.test2 = Node('test2')
        self.cons = Connections(None, 'cons')
        self.cons += self.test1
        self.cons += self.test2

    def test_01_properties(self):
        self.assertListEqual(sorted(self.cons.names),
                             ['test1', 'test2'])
        self.assertListEqual(sorted(self.cons.slaves),
                             [self.test1, self.test2])
    def test_02_contains(self):
        self.assertTrue('test1' in self.cons)
        self.assertTrue('test2' in self.cons)
        self.assertTrue(self.test1 in self.cons)
        self.assertTrue(self.test2 in self.cons)
    def test_03_iterable(self):
        devices = []
        for device in self.cons:
            devices.append(device)
        self.assertListEqual(devices, [self.test1, self.test2])

class Test01Self2Node(unittest.TestCase):

    def tearDown(self):
        Element.clear_registry()

    def test_01_iadd(self):
        test = Connections(None, 'test')
        n1 = Node('n1', 'Q')
        test += n1
        self.assertIsInstance(test, Connections)
        self.assertIsInstance(test.n1, Node)
        self.assertIs(test.n1, n1)
        test += n1
        self.assertIsInstance(test, Connections)
        self.assertIsInstance(test.n1, Node)
        self.assertIs(test.n1, n1)
        n2 = Node('n2', 'T')
        test += n2
        self.assertIs(test.n1, n1)
        self.assertIs(test.n2, n2)

    def test_02_variables(self):
        test = Connections(None, 'test')
        self.assertListEqual(test.variables, [])
        n1 = Node('n1', 'Q')
        test += n1
        self.assertListEqual(test.variables, ['Q'])
        n2 = Node('n2', 'T')
        test += n2
        self.assertListEqual(sorted(test.variables), ['Q', 'T'])
        with self.assertRaises(AttributeError):
            test.X

# import...
# ...from standard library
from __future__ import division, print_function
import unittest
# ...from HydPy
from hydpy.core.devicetools import *
from hydpy.core.selectiontools import *

class Test01SelectionInitialization(unittest.TestCase):

    def tearDown(self):
        Node.clear_registry()
        Element.clear_registry()

    def test_01_empty(self):
        test = Selection('test')
        self.assertIsInstance(test, Selection)
        self.assertIsInstance(test.nodes, Nodes)
        self.assertIsInstance(test.elements, Elements)
        self.assertFalse(bool(test))
    def test_02_nonempty(self):
        test = Selection('test', nodes='a', elements=('b', 'c'))
        self.assertIsInstance(test, Selection)
        self.assertIsInstance(test.nodes, Nodes)
        self.assertIsInstance(test.elements, Elements)
        self.assertEqual(test.nodes.a.name, 'a')
        self.assertEqual(test.elements.b.name, 'b')
        self.assertEqual(test.elements.c.name, 'c')
        self.assertTrue(bool(test))
    def test_03_copy(self):
        test1 = Selection('test1', nodes='a', elements=('b', 'c'))
        test2 = test1.copy('test2')
        self.assertIsInstance(test2, Selection)
        self.assertEqual(test2.name, 'test2')
        self.assertIsNot(test1, test2)
        self.assertEqual(test2.nodes, test1.nodes)
        self.assertEqual(test2.elements, test1.elements)


class Test02SelectionSelect(unittest.TestCase):

    def setUp(self):
        # e1 + e2 -> n_Q1 -> e3 -> n_Q2 -> e4 + (e7 -> n_Q4 -> e6)-> n_Q3
        # e1 + e2 -> n_T1 -> e3 -> n_T2 -> e5 -> n_T3
        self.n_Q1 = Node('n_Q1', 'Q')
        self.n_T1 = Node('n_T1', 'T')
        self.e1 = Element('e1', outlets='n_Q1')
        self.e1 = Element('e1', outlets='n_T1')
        self.e2 = Element('e2', outlets='n_Q1')
        self.e2 = Element('e2', outlets='n_T1')
        self.n_Q2 = Node('n_Q2', 'Q')
        self.n_T2 = Node('n_T2', 'T')
        self.e3 = Element('e3', inlets='n_Q1', outlets='n_Q2')
        self.e3 = Element('e3', inlets='n_T1', outlets='n_T2')
        self.n_Q3 = Node('n_Q3', 'Q')
        self.e4 = Element('e4', inlets='n_Q2', outlets='n_Q3')
        self.n_T3 = Node('n_T3', 'T')
        self.e5 = Element('e5', inlets='n_T2', outlets='n_T3')
        self.n_Q4 = Node('n_Q4', 'Q')
        self.e6 = Element('e6', inlets='n_Q4', outlets='n_Q3')
        self.e7 = Element('e7', outlets='n_Q4')

        self.complete = Selection('complete',
                                  nodes=Node.registered_nodes(),
                                  elements=Element.registered_elements())

    def tearDown(self):
        Node.clear_registry()
        Element.clear_registry()

    def test_01_nextelement(self):
        nodes, elements = self.complete._nextelement(self.e1,
                                                  Nodes(), Elements())
        self.assertEqual(nodes, Nodes())
        self.assertEqual(elements, Elements(self.e1))
    def test_02_nextnode(self):
        nodes, elements = self.complete._nextnode(self.n_Q1,
                                                  Nodes(), Elements())
        self.assertEqual(nodes, Nodes(self.n_Q1))
        self.assertEqual(elements, Elements(self.e1, self.e2))
    def test_03_select_upstream(self):
        test = self.complete.copy('test').select_upstream(self.n_Q3)
        reference = self.complete.copy('reference')
        del(reference.elements.e5)
        del(reference.nodes.n_T2)
        del(reference.nodes.n_T3)
        self.assertEqual(test.nodes, reference.nodes)
        self.assertEqual(test.elements, reference.elements)
    def test_04_deselect_upstream(self):
        test = self.complete.copy('test').deselect_upstream(self.n_Q3)
        reference = Selection('reference', ['n_T2', 'n_T3'], 'e5')
        self.assertEqual(test.nodes, reference.nodes)
        self.assertEqual(test.elements, reference.elements)

#    def test_05_selectmodelclasses(self):
#        with self.assertRaises(RuntimeError):
#            self.complete.copy('test').getby_modelclasses('HBV96_zone')

    def test_06_select_nodenames(self):
        test = self.complete.copy('test').select_nodenames('n_Q1', 'n_T', 'NO')
        reference = self.complete.copy('test')
        reference.nodes = Nodes('n_Q1', 'n_T1', 'n_T2', 'n_T3')
        self.assertEqual(test.nodes, reference.nodes)
        self.assertEqual(test.elements, reference.elements)
    def test_07_deselect_nodenames(self):
        test = self.complete.copy('test').deselect_nodenames('n_Q1', 'n_T',
                                                             'NO')
        reference = self.complete.copy('test')
        del(reference.nodes.n_Q1)
        del(reference.nodes.n_T1)
        del(reference.nodes.n_T2)
        del(reference.nodes.n_T3)
        self.assertEqual(test.nodes, reference.nodes)
        self.assertEqual(test.elements, reference.elements)
    def test_08_select_elementnames(self):
        test = self.complete.copy('test').select_elementnames('e')
        self.assertEqual(test.nodes, self.complete.nodes)
        self.assertEqual(test.elements, self.complete.elements)
        test = self.complete.copy('test').select_elementnames('2')
        reference = self.complete.copy('test')
        reference.elements = Elements('e2')
    def test_09_deselect_elementnames(self):
        test = self.complete.copy('test').deselect_elementnames('e')
        self.assertEqual(test.nodes, self.complete.nodes)
        self.assertEqual(test.elements, Elements())
        test = self.complete.copy('test').deselect_elementnames('2')
        reference = self.complete.copy('test')
        del(reference.elements.e2)
        self.assertEqual(test.nodes, reference.nodes)
        self.assertEqual(test.elements, reference.elements)


class Test03SelectionMagic(unittest.TestCase):

    def tearDown(self):
        Node.clear_registry()
        Element.clear_registry()

    def test_01_len(self):
        test1 = Selection('test1', ['n1', 'n2'], [])
        self.assertEqual(len(test1), 2)
        test2 = Selection('test2', [], ['e1'])
        self.assertEqual(len(test2), 1)

    def test_02_iadd(self):
        test = Selection('sel',
                         ['n1', 'n2', 'n3'],
                         ['e1', 'e2', 'e3'])
        test += Selection('add', ['n3', 'n4'], [])
        reference = Selection('sel',
                              ['n1', 'n2', 'n3', 'n4'],
                              ['e1', 'e2', 'e3'])
        self.assertEqual(test.nodes, reference.nodes)
        self.assertEqual(test.elements, reference.elements)
        test += Selection('add', [], ['e4', 'e5'])
        reference = Selection('sel',
                              ['n1', 'n2', 'n3', 'n4'],
                              ['e1', 'e2', 'e3', 'e4', 'e5'])
        self.assertEqual(test.nodes, reference.nodes)
        self.assertEqual(test.elements, reference.elements)
    def test_03_isub(self):
        test = Selection('sel',
                         ['n1', 'n2', 'n3'],
                         ['e1', 'e2', 'e3'])
        test -= Selection('sub', ['n3', 'n4'], [])
        reference = Selection('sel',
                              ['n1', 'n2'],
                              ['e1', 'e2', 'e3'])
        test -= Selection('sub', [], ['e1', 'e2', 'e3'])
        reference = Selection('sel',
                              ['n1', 'n2'],
                              [])
        self.assertEqual(test.nodes, reference.nodes)
        self.assertEqual(test.elements, reference.elements)

# -*- coding: utf-8 -*-
"""
Created on Thu Dec 01 10:53:44 2016

@author: tyralla
"""

# import...
# ...from standard library
from __future__ import division, print_function
import unittest
# ...from site-packages
import numpy
# ...from HydPy
from hydpy import pub
from hydpy.core import indextools
from hydpy.core import timetools

class Test01MonthOfYear(unittest.TestCase):

    def setUp(self):
        self.indexer = indextools.Indexer()

    def tearDown(self):
        pub.timegrids = None

    def test_01_manual_mode(self):
        with self.assertRaises(BaseException):
            self.indexer.monthofyear = 'a'
        with self.assertRaises(BaseException):
            self.indexer.monthofyear = ['a', 'b']
        with self.assertRaises(ValueError):
            self.indexer.monthofyear = [[1, 2], [3, 4]]
        self.indexer.monthofyear = [1,2]
        self.assertIsInstance(self.indexer.monthofyear, numpy.ndarray)
        self.assertTupleEqual(tuple(self.indexer.monthofyear), (1, 2))
        del(self.indexer.monthofyear)
        self.assertIsNone(self.indexer._monthofyear)
        pub.timegrids = timetools.Timegrids(timetools.Timegrid('01.01.2004',
                                                               '1.01.2005',
                                                               '1d'))
        with self.assertRaises(ValueError):
            self.indexer.monthofyear = [1,2]

    def test_02_automatic_mode(self):
        with self.assertRaises(RuntimeError):
            self.indexer.monthofyear
        pub.timegrids = timetools.Timegrids(timetools.Timegrid('01.01.2004',
                                                               '1.01.2005',
                                                               '1d'))
        self.assertIsInstance(self.indexer.monthofyear, numpy.ndarray)
        self.assertEqual(len(self.indexer.monthofyear), 366)
        self.assertTupleEqual(tuple(self.indexer.monthofyear),
                              tuple(31*[0]+29*[1]+31*[2]+30*[3]+
                                    31*[4]+30*[5]+31*[6]+31*[7]+
                                    30*[8]+31*[9]+30*[10]+31*[11]))
        self.assertIs(self.indexer.monthofyear, self.indexer.monthofyear)


class Test02DayOfYear(unittest.TestCase):

    def setUp(self):
        self.indexer = indextools.Indexer()

    def tearDown(self):
        pub.timegrids = None

    def test_01_manual_mode(self):
        with self.assertRaises(BaseException):
            self.indexer.dayofyear = 'a'
        with self.assertRaises(BaseException):
            self.indexer.dayofyear = ['a', 'b']
        with self.assertRaises(ValueError):
            self.indexer.dayofyear = [[1, 2], [3, 4]]
        self.indexer.dayofyear = [1,2]
        self.assertIsInstance(self.indexer.dayofyear, numpy.ndarray)
        self.assertTupleEqual(tuple(self.indexer.dayofyear), (1, 2))
        del self.indexer.dayofyear
        self.assertIsNone(self.indexer._dayofyear)
        pub.timegrids = timetools.Timegrids(timetools.Timegrid('01.01.2004',
                                                               '1.01.2005',
                                                               '1d'))
        with self.assertRaises(ValueError):
            self.indexer.dayofyear = [1,2]

    def test_02_automatic_mode(self):
        with self.assertRaises(RuntimeError):
            self.indexer.dayofyear
        pub.timegrids = timetools.Timegrids(timetools.Timegrid('01.01.2004',
                                                               '1.01.2005',
                                                               '1d'))
        self.assertIsInstance(self.indexer.dayofyear, numpy.ndarray)
        self.assertEqual(len(self.indexer.dayofyear), 366)
        self.assertTupleEqual(tuple(self.indexer.dayofyear),
                              tuple(range(366)))
        pub.timegrids = timetools.Timegrids(timetools.Timegrid('01.01.2005',
                                                               '1.01.2006',
                                                               '1d'))
        del self.indexer.dayofyear
        self.assertIsInstance(self.indexer.dayofyear, numpy.ndarray)
        self.assertEqual(len(self.indexer.dayofyear), 365)
        self.assertTupleEqual(
           tuple(self.indexer.dayofyear),
           tuple(list(range(31+28))+list(range(31+28+1, 366)))
           )
        self.assertIs(self.indexer.dayofyear, self.indexer.dayofyear)

class Test03TimeOfYear(unittest.TestCase):

    def setUp(self):
        self.indexer = indextools.Indexer()

    def tearDown(self):
        pub.timegrids = None

    def test_01_manual_mode(self):
        with self.assertRaises(BaseException):
            self.indexer.timeofyear = 'a'
        with self.assertRaises(BaseException):
            self.indexer.timeofyear = ['a', 'b']
        with self.assertRaises(ValueError):
            self.indexer.timeofyear = [[1, 2], [3, 4]]
        self.indexer.timeofyear = [1, 2]
        self.assertIsInstance(self.indexer.timeofyear, numpy.ndarray)
        self.assertTupleEqual(tuple(self.indexer.timeofyear), (1, 2))
        del self.indexer.timeofyear
        self.assertIsNone(self.indexer._timeofyear)
        pub.timegrids = timetools.Timegrids(timetools.Timegrid('01.01.2004',
                                                               '1.01.2005',
                                                               '1d'))
        with self.assertRaises(ValueError):
            self.indexer.timeofyear = [1, 2]

    def test_02_automatic_mode(self):
        with self.assertRaises(RuntimeError):
            self.indexer.timeofyear
        pub.timegrids = timetools.Timegrids(timetools.Timegrid('01.01.2004',
                                                               '1.01.2005',
                                                               '1d'))
        self.assertIsInstance(self.indexer.timeofyear, numpy.ndarray)
        self.assertEqual(len(self.indexer.timeofyear), 366)
        self.assertTupleEqual(tuple(self.indexer.timeofyear),
                              tuple(range(366)))
        pub.timegrids = timetools.Timegrids(timetools.Timegrid('01.01.2005',
                                                               '1.01.2006',
                                                               '1d'))
        del self.indexer.timeofyear
        self.assertIsInstance(self.indexer.timeofyear, numpy.ndarray)
        self.assertEqual(len(self.indexer.timeofyear), 365)
        self.assertTupleEqual(
           tuple(self.indexer.timeofyear),
           tuple(list(range(31+28))+list(range(31+28+1, 366)))
           )
        self.assertIs(self.indexer.timeofyear, self.indexer.timeofyear)

"""
HydPy

An interactive framework for the developement and a application of
hydrological models.
"""
# import...
# ...from standard library
import sys
import warnings
try:
    import builtins
except ImportError:
    import __builtin__ as builtins
# ...from site-packages
from hydpy.core.exceptiontools import OptionalImport
import numpy
from numpy import nan
from scipy import integrate
pandas = OptionalImport('import pandas')
pyplot = OptionalImport('from matplotlib import pyplot')
# ...from HydPy
from hydpy import pub
from hydpy.core import dummytools
from hydpy.core import indextools
from hydpy.core import optiontools
from hydpy.cythons import configutils
from hydpy.core.autodoctools import prepare_mainsubstituter
from hydpy.core.auxfiletools import Auxfiler
from hydpy.core.devicetools import Element
from hydpy.core.devicetools import Elements
from hydpy.core.devicetools import Node
from hydpy.core.devicetools import Nodes
from hydpy.core.hydpytools import HydPy
from hydpy.core.importtools import prepare_model
from hydpy.core.importtools import reverse_model_wildcard_import
from hydpy.core.objecttools import HydPyDeprecationWarning
from hydpy.core.objecttools import print_values
from hydpy.core.objecttools import round_
from hydpy.core.selectiontools import Selection
from hydpy.core.selectiontools import Selections
from hydpy.core.timetools import Date
from hydpy.core.timetools import Period
from hydpy.core.timetools import Timegrid
from hydpy.core.timetools import Timegrids
from hydpy.core.testtools import IntegrationTest
from hydpy.core.testtools import Open
from hydpy.core.testtools import UnitTest
from hydpy.auxs.armatools import ARMA
from hydpy.auxs.armatools import MA
from hydpy.auxs.anntools import ANN
from hydpy.auxs.anntools import ann
from hydpy.auxs.anntools import SeasonalANN
from hydpy.auxs.iuhtools import LinearStorageCascade
from hydpy.auxs.iuhtools import TranslationDiffusionEquation
from hydpy.auxs.networktools import RiverBasinNumber
from hydpy.auxs.networktools import RiverBasinNumbers
from hydpy.auxs.networktools import RiverBasinNumbers2Selection
from hydpy.auxs.statstools import bias_abs
from hydpy.auxs.statstools import bias_rel
from hydpy.auxs.statstools import calc_mean_time
from hydpy.auxs.statstools import calc_mean_time_deviation
from hydpy.auxs.statstools import corr
from hydpy.auxs.statstools import evaluationtable
from hydpy.auxs.statstools import hsepd
from hydpy.auxs.statstools import hsepd_manual
from hydpy.auxs.statstools import hsepd_pdf
from hydpy.auxs.statstools import nse
from hydpy.auxs.statstools import prepare_arrays
from hydpy.auxs.statstools import std_ratio


pub.options = optiontools.Options()
pub.indexer = indextools.Indexer()
pub.config = configutils.Config()
dummies = dummytools.Dummies()   # pylint: disable=invalid-name


def customwarn(message, category, filename, lineno, file=None, line=None):
    """Redirect warnings to `stdout`."""
    sys.stdout.write(warnings.formatwarning(
        message, category, filename, lineno))


warnings.showwarning = customwarn
warnings.filterwarnings('always', category=HydPyDeprecationWarning)
warnings.filterwarnings('ignore', r'All-NaN (slice|axis) encountered')
warnings.filterwarnings('error', category=integrate.IntegrationWarning)

# Numpy introduced new string representations in version 1.14 affecting
# our doctests.  Hence, the old style is selected for now:
try:
    # pylint: disable=unexpected-keyword-arg
    numpy.set_printoptions(legacy='1.13')
except TypeError:
    pass

substituter = prepare_mainsubstituter()


__all__ = ['builtins',
           'pandas',
           'pyplot',
           'pub',
           'Auxfiler',
           'Element',
           'Elements',
           'Node',
           'Nodes',
           'HydPy',
           'prepare_model',
           'reverse_model_wildcard_import',
           'print_values',
           'round_',
           'Selection',
           'Selections',
           'Date',
           'Period',
           'Timegrid',
           'Timegrids',
           'IntegrationTest',
           'Open',
           'UnitTest',
           'ARMA',
           'MA',
           'ANN',
           'ann',
           'SeasonalANN',
           'LinearStorageCascade',
           'TranslationDiffusionEquation',
           'RiverBasinNumber',
           'RiverBasinNumbers',
           'RiverBasinNumbers2Selection',
           'nan',
           'bias_abs',
           'bias_rel',
           'calc_mean_time',
           'calc_mean_time_deviation',
           'corr',
           'evaluationtable',
           'hsepd',
           'hsepd_manual',
           'hsepd_pdf',
           'nse',
           'prepare_arrays',
           'std_ratio']

# -*- coding: utf-8 -*-
"""Module |config| allows the user to configure `HydPy`.

The available options should not be changed during runtime.
"""
use_autodoc = True
"""Flag that indicates whether HydPy's automatic documentation manipulation
features should be applyied or not.  It is imperative to set it to `True`
before one uses `Sphinx` to generate the online documentation.  However,
one can set it to `False` for regular applications, which reduces HydPy's
initialisation time."""

# import from...
# ...the standard library
import sys as _sys

projectname = None

options = None
indexer = None
networkmanager = None
controlmanager = None
conditionmanager = None
sequencemanager = None
timegrids = None

pyversion = int(_sys.version[0])
_printprogress_indentation = -4

_is_hydpy_bundled = getattr(_sys, 'frozen', False)
"""This parameter is set `True` within HydPy executables only.
Then different features that do not make much sense within an
executable (e.g. the cythonization of models) are disabled.
"""

#

# -*- coding: utf-8 -*-
"""This module implements rudimantary artificial neural network tools,
required for some models implemented in the HydPy framework.

A note for developers: some of the implemented features are to be applied
during model simulations are in some other way performance-critical.  Hence,
the actual calculations are defined in the Cython extension module
|annutils|.
"""

# import...
# ...from standard library
from __future__ import division, print_function
# ...from site-packages
import numpy
from hydpy import pyplot
# ...from HydPy
from hydpy.core import abctools
from hydpy.core import autodoctools
from hydpy.core import exceptiontools
from hydpy.core import objecttools
from hydpy.core import parametertools
from hydpy.core import timetools
from hydpy.cythons import annutils   # pylint: disable=no-name-in-module


class ANN(object):
    """Multi-layer feed forward artificial neural network.

    The applied activation function is the logistic function:

      :math:`f(x) = \\frac{1}{1+exp(-x)}`

    Class |anntools.ANN| is intended to be subclassed for the derivation of
    very complex control parameters.  Its original purpose was to allow for
    defining arbitrary continuous relationsships between the water stored
    in a dam and the associated water stage (see model ...).  However,
    class |anntools.ANN| can also be applied directly, as shown in the
    following examples.  But if you are looking for a flexible stand-alone
    artifical neural network implementation in Python, you will find much
    more general tools easily.

    Firstly, define the most single artificial neural network consisting of
    only one input node, neuron, and output node respectively, and pass
    some arbitrary network parameters:

    >>> from hydpy import ANN, nan
    >>> ann = ANN()
    >>> ann(nmb_inputs=1, nmb_neurons=(1,), nmb_outputs=1,
    ...     weights_input=4.0, weights_output=3.0,
    ...     intercepts_hidden=-16.0, intercepts_output=-1.0)

    The following loop subsequently sets the values 0 to 8 as input values,
    performs the calculateion, and prints out the final output.  As to be
    expected, the results show the shape of the logistic function:

    >>> from hydpy import round_
    >>> for input_ in range(9):
    ...     ann.inputs[0] = input_
    ...     ann.process_actual_input()
    ...     round_([input_, ann.outputs[0]])
    0, -1.0
    1, -0.999982
    2, -0.998994
    3, -0.946041
    4, 0.5
    5, 1.946041
    6, 1.998994
    7, 1.999982
    8, 2.0

    One can also directly plot the resulting graph:

    >>> ann.plot(0.0, 8.0)

    .. testsetup::

        >>> from matplotlib import pyplot
        >>> pyplot.close()

    The following example shows that everything works well for more complex
    single layer networks also (manual tests have been performed in a
    spreadsheet program):

    >>> ann.nmb_inputs = 3
    >>> ann.nmb_neurons = (4,)
    >>> ann.nmb_outputs = 2
    >>> ann.weights_input = [[ 0.2, -0.1, -1.7,  0.6],
    ...                      [ 0.9,  0.2,  0.8,  0.0],
    ...                      [-0.5, -1.0,  2.3, -0.4]]
    >>> ann.weights_output = [[ 0.0,  2.0],
    ...                       [-0.5,  1.0],
    ...                       [ 0.4,  2.4],
    ...                       [ 0.8, -0.9]]
    >>> ann.intercepts_hidden = [ 0.9,  0.0, -0.4, -0.2]
    >>> ann.intercepts_output = [ 1.3, -2.0]
    >>> ann.inputs = [-0.1,  1.3,  1.6]
    >>> ann.process_actual_input()
    >>> round_(ann.outputs)
    1.822222, 1.876983

    The next example shows how to solve the XOR problem with a two layer
    network.  As usual, `1` stands for `True` and `0` stands for `False`.

    We define a network with two inputs (`I1` and `I2`), two neurons in
    the first hidden layer (`H11` and `H12`), one neuron in the second
    hidden layer (`H2`), and a single output (`O1`):

    >>> ann.nmb_inputs = 2
    >>> ann.nmb_neurons = (2, 1)
    >>> ann.nmb_outputs = 1

    The value of `O1` shall be identical with the activation of `H2`:

    >>> ann.weights_output = 1.0
    >>> ann.intercepts_output = 0.0

    All intercepts of the neurons of the hidden layer are set to 750,
    so that an input of 500 results in an activation of approximately
    zero and an input of 1000 results in an activation of approximately
    one (note that matrix entries are not required should preferably be
    initialized with `nan` to avoid confusion):

    >>> ann.intercepts_hidden = [[-750.0, -750.0],
    ...                          [-750.0, nan]]

    The weighting factor between the both inputs and `H11` is 1000.
    Hence, one `True` input is sufficient to activate `H1`.  In contrast,
    the weighting factor between the both inputs and `H12` is 500 only.
    Hence, two `True` inputs are required to activate `H12`:

    >>> ann.weights_input= [[1000.0, 500.0],
    ...                     [1000.0, 500.0]]

    The weighting factor between `H11` and `H2` is 1000.  Hence, in
    principle, `H11` can activate `H2`.  However, the weighting factor
    between `H12` and `H2` is -1000.  Hence, `H12` is able to prevent
    `H2` from becoming activated even when `H11` is activated:

    >>> ann.weights_hidden= [[[1000.0, nan],
    ...                      [-1000.0, nan]]]

    To recapitulate, `H11` determines if at least one input is `True`,
    `H12` determines if both inputs are `True`, and `H2` determines
    if exactly one input is `True`, which is the solution for the XOR-problem:

    >>> ann
    ann(nmb_inputs=2,
        nmb_neurons=(2, 1),
        nmb_outputs=1,
        weights_input=[[1000.0, 500.0],
                       [1000.0, 500.0]],
        weights_hidden=[[[1000.0, nan],
                         [-1000.0, nan]]],
        weights_output=[[1.0]],
        intercepts_hidden=[[-750.0, -750.0],
                           [-750.0, nan]],
        intercepts_output=[0.0])

    The following calculation confirms that the network is properly
    configured:

    >>> for inputs in ((0.0, 0.0),
    ...                (1.0, 0.0),
    ...                (0.0, 1.0),
    ...                (1.0, 1.0)):
    ...    ann.inputs = inputs
    ...    ann.process_actual_input()
    ...    print(inputs[0], inputs[1], ann.outputs[0])
    0.0 0.0 0.0
    1.0 0.0 1.0
    0.0 1.0 1.0
    1.0 1.0 0.0

    To elaborate on the last calculation, the corresponding activations
    of the hidden neurons are shown. As both inputs are `True`, both
    `H12` (upper left value) and `H22` (upper right value) activated,
    but `H2` (lower left value) is not:

    >>> ann.neurons
    array([[ 1.,  1.],
           [ 0.,  0.]])

    The last defined configuration is used in some examples of the
    documentation of the members of class |anntools.ANN|:

    >>> from hydpy import dummies
    >>> dummies.ann = ann

    Note that Python class |anntools.ANN| handles a corresponding
    Cython extension class defined in |annutils|, which does not
    protect itself against segmentation faults. But class
    |anntools.ANN| takes up this task, meaning using its public
    members should always result in readable exceptions instead of
    program crashes, e.g.:

    >>> ANN().nmb_layers
    Traceback (most recent call last):
    ...
    AttributeNotReady: Attribute `nmb_layers` of object `ann` \
is not usable so far.
    """
    NDIM = 0
    TYPE = 'annutils.ANN'
    TIME = None
    SPAN = (None, None)

    parameterstep = parametertools.Parameter.__dict__['parameterstep']
    simulationstep = parametertools.Parameter.__dict__['simulationstep']

    def __init__(self):
        self.subpars = None
        self.fastaccess = objecttools.FastAccess()
        self._isready = exceptiontools.IsReady(
            false=('nmb_inputs',  'nmb_outputs', 'nmb_neurons'))
        self._cann = annutils.ANN()
        self._max_nmb_neurons = None

    def connect(self, subpars):
        """Connect the actual |anntools.ANN| object with the given
        |SubParameters| object."""
        self.subpars = subpars
        self.fastaccess = subpars.fastaccess
        setattr(self.fastaccess, self.name, self._cann)

    name = property(objecttools.name)

    def __call__(self, nmb_inputs=1, nmb_neurons=(1,), nmb_outputs=1,
                 weights_input=None, weights_output=None, weights_hidden=None,
                 intercepts_hidden=None, intercepts_output=None):
        self.nmb_inputs = nmb_inputs
        self.nmb_outputs = nmb_outputs
        self.nmb_neurons = nmb_neurons
        self.weights_input = weights_input
        self.weights_hidden = weights_hidden
        self.weights_output = weights_output
        self.intercepts_hidden = intercepts_hidden
        self.intercepts_output = intercepts_output
        del self.inputs
        del self.outputs
        del self.neurons

    def _update_shapes(self):
        if self._isready:
            del self.weights_input
            del self.weights_hidden
            del self.weights_output
            del self.intercepts_hidden
            del self.intercepts_output
            del self.inputs
            del self.outputs
            del self.neurons

    def _get_nmb_inputs(self):
        """Number of input nodes.

        >>> from hydpy import dummies
        >>> dummies.ann.nmb_inputs
        2
        """
        return self._cann.nmb_inputs

    def _set_nmb_inputs(self, value):
        self._cann.nmb_inputs = int(value)
        self._update_shapes()

    nmb_inputs = exceptiontools.protected_property(
        'nmb_inputs', _get_nmb_inputs, _set_nmb_inputs)

    def _get_nmb_outputs(self):
        """Number of output nodes.

        >>> from hydpy import dummies
        >>> dummies.ann.nmb_outputs
        1
        """
        return self._cann.nmb_outputs

    def _set_nmb_outputs(self, value):
        self._cann.nmb_outputs = int(value)
        self._update_shapes()

    nmb_outputs = exceptiontools.protected_property(
        'nmb_outputs', _get_nmb_outputs, _set_nmb_outputs)

    def _get_nmb_layers(self):
        """Number of hidden layers.

        >>> from hydpy import dummies
        >>> dummies.ann.nmb_layers
        2
        """
        return self._cann.nmb_layers

    nmb_layers = exceptiontools.dependent_property(
        'nmb_layers', _get_nmb_layers)

    def _get_nmb_neurons(self):
        """Number of neurons of the hidden layers.

        >>> from hydpy import dummies
        >>> dummies.ann.nmb_neurons
        (2, 1)
        """
        return tuple(numpy.asarray(self._cann.nmb_neurons))

    def _set_nmb_neurons(self, value):
        self._cann.nmb_neurons = numpy.array(value, dtype=int, ndmin=1)
        self._cann.nmb_layers = len(value)
        self._max_nmb_neurons = max(value)
        self._update_shapes()

    nmb_neurons = exceptiontools.protected_property(
        'nmb_neurons', _get_nmb_neurons, _set_nmb_neurons)

    def _get_weights_input(self):
        """Weights between all input nodes and neurons of the first hidden
        layer.

        The input nodes and the neurons are varied on the first axis and
        on the second axis of the 2-dimensional array:

        >>> from hydpy import ANN
        >>> ann = ANN()
        >>> ann(nmb_inputs=2, nmb_neurons=(3,))
        >>> ann.weights_input
        array([[ 0.,  0.,  0.],
               [ 0.,  0.,  0.]])

        It is allowed to set values via slicing:

        >>> ann.weights_input[:, 0] = 1.
        >>> ann.weights_input
        array([[ 1.,  0.,  0.],
               [ 1.,  0.,  0.]])

        If possible, type conversions are performed:

        >>> ann.weights_input = '2'
        >>> ann.weights_input
        array([[ 2.,  2.,  2.],
               [ 2.,  2.,  2.]])

        One can assign whole matrices directly:

        >>> import numpy
        >>> ann.weights_input = numpy.eye(2, 3)
        >>> ann.weights_input
        array([[ 1.,  0.,  0.],
        ...    [ 0.,  1.,  0.]])

        One can also delete the values contained in the array:

        >>> del ann.weights_input
        >>> ann.weights_input
        array([[ 0.,  0.,  0.],
        ...    [ 0.,  0.,  0.]])

        Errors like wrong shapes (or unconvertible inputs) result in error
        messages:

        >>> ann.weights_input = numpy.eye(3)
        Traceback (most recent call last):
        ...
        ValueError: While trying to set the input weights of the artificial \
neural network `ann` of element `?`, the following error occured: could not \
broadcast input array from shape (3,3) into shape (2,3)
        """
        return numpy.asarray(self._cann.weights_input)

    def _set_weights_input(self, values):
        if values is None:
            self._del_weights_input()
        else:
            try:
                self._cann.weights_input = numpy.full(self.shape_weights_input,
                                                      values, dtype=float)
            except BaseException:
                objecttools.augment_excmessage(
                    'While trying to set the input weights of the artificial '
                    'neural network `%s` of element `%s`'
                    % (self.name, objecttools.devicename(self)))

    def _del_weights_input(self):
        self._cann.weights_input = numpy.zeros(self.shape_weights_input)

    weights_input = exceptiontools.dependent_property(
        'weights_input', _get_weights_input,
        _set_weights_input, _del_weights_input)

    @property
    def shape_weights_input(self):
        """Shape of the array containing the input weights.

        >>> from hydpy import dummies
        >>> dummies.ann.shape_weights_input
        (2, 2)
        """
        return (self.nmb_inputs, self.nmb_neurons[0])

    @property
    def nmb_weights_input(self):
        """Number of input weights.

        >>> from hydpy import dummies
        >>> dummies.ann.nmb_weights_input
        4
        """
        return self.nmb_neurons[0]*self.nmb_inputs

    def _get_weights_output(self):
        """Weights between all neurons of the last hidden layer and the output
        nodes.

        The neurons and the output nodes are varied on the first axis and
        on the second axis of the 2-dimensional array:

        >>> from hydpy import ANN
        >>> ann = ANN()
        >>> ann(nmb_outputs=2, nmb_neurons=(3,))
        >>> ann.weights_output
        array([[ 0.,  0.],
               [ 0.,  0.],
               [ 0.,  0.]])

        It is allowed to set values via slicing:

        >>> ann.weights_output[:, 0] = 1.
        >>> ann.weights_output
        array([[ 1.,  0.],
               [ 1.,  0.],
               [ 1.,  0.]])

        If possible, type conversions are performed:

        >>> ann.weights_output = '2'
        >>> ann.weights_output
        array([[ 2.,  2.],
               [ 2.,  2.],
               [ 2.,  2.]])

        One can assign whole matrices directly:

        >>> import numpy
        >>> ann.weights_output = numpy.eye(3, 2)
        >>> ann.weights_output
        array([[ 1.,  0.],
               [ 0.,  1.],
               [ 0.,  0.]])

        One can also delete the values contained in the array:

        >>> del ann.weights_output
        >>> ann.weights_output
        array([[ 0.,  0.],
               [ 0.,  0.],
               [ 0.,  0.]])

        Errors like wrong shapes (or unconvertible inputs) result in error
        messages:

        >>> ann.weights_output = numpy.eye(3)
        Traceback (most recent call last):
        ...
        ValueError: While trying to set the output weights of the artificial \
neural network `ann` of element `?`, the following error occured: could not \
broadcast input array from shape (3,3) into shape (3,2)
        """
        return numpy.asarray(self._cann.weights_output)

    def _set_weights_output(self, values):
        if values is None:
            self._del_weights_output()
        else:
            try:
                self._cann.weights_output = numpy.full(
                    self.shape_weights_output, values, dtype=float)
            except BaseException:
                objecttools.augment_excmessage(
                    'While trying to set the output weights of the artificial '
                    'neural network `%s` of element `%s`'
                    % (self.name, objecttools.devicename(self)))

    def _del_weights_output(self):
        self._cann.weights_output = numpy.zeros(self.shape_weights_output)

    weights_output = exceptiontools.dependent_property(
        'weights_output',  _get_weights_output,
        _set_weights_output, _del_weights_output)

    @property
    def shape_weights_output(self):
        """Shape of the array containing the output weights.

        >>> from hydpy import dummies
        >>> dummies.ann.shape_weights_output
        (1, 1)
        """
        return (self.nmb_neurons[-1], self.nmb_outputs)

    @property
    def nmb_weights_output(self):
        """Number of output weights.

        >>> from hydpy import dummies
        >>> dummies.ann.nmb_weights_output
        1
        """
        return self.nmb_neurons[-1]*self.nmb_outputs

    def _get_weights_hidden(self):
        """Weights between between the neurons of the different hidden layers.

        The layers are varied on the first axis, the neurons of the respective
        upstream layer on the second axis and the neurons of the respective
        downstream layer on the third axis of a 3-dimensional array:

        >>> from hydpy import ANN
        >>> ann = ANN()
        >>> ann(nmb_neurons=(3, 2, 3))
        >>> ann.weights_hidden
        array([[[  0.,   0.,  nan],
                [  0.,   0.,  nan],
                [  0.,   0.,  nan]],
        <BLANKLINE>
               [[  0.,   0.,   0.],
                [  0.,   0.,   0.],
                [ nan,  nan,  nan]]])

        It is allowed to set values via slicing:

        >>> ann.weights_hidden[1, :, 0] = 1.
        >>> ann.weights_hidden
        array([[[  0.,   0.,  nan],
                [  0.,   0.,  nan],
                [  0.,   0.,  nan]],
        <BLANKLINE>
               [[  1.,   0.,   0.],
                [  1.,   0.,   0.],
                [  1.,  nan,  nan]]])

        If possible, type conversions are performed:

        >>> ann.weights_hidden = '2'
        >>> ann.weights_hidden
        array([[[ 2.,  2.,  2.],
                [ 2.,  2.,  2.],
                [ 2.,  2.,  2.]],
        <BLANKLINE>
               [[ 2.,  2.,  2.],
                [ 2.,  2.,  2.],
                [ 2.,  2.,  2.]]])

        One can assign whole matrices directly:

        >>> import numpy
        >>> ann.weights_hidden = numpy.eye(3)
        >>> ann.weights_hidden
        array([[[ 1.,  0.,  0.],
                [ 0.,  1.,  0.],
                [ 0.,  0.,  1.]],
        <BLANKLINE>
               [[ 1.,  0.,  0.],
                [ 0.,  1.,  0.],
                [ 0.,  0.,  1.]]])

        One can also delete the values contained in the array:

        >>> del ann.weights_hidden
        >>> ann.weights_hidden
        array([[[  0.,   0.,  nan],
                [  0.,   0.,  nan],
                [  0.,   0.,  nan]],
        <BLANKLINE>
               [[  0.,   0.,   0.],
                [  0.,   0.,   0.],
                [ nan,  nan,  nan]]])

        Errors like wrong shapes (or unconvertible inputs) result in error
        messages:

        >>> ann.weights_hidden = numpy.eye(3, 2)
        Traceback (most recent call last):
        ...
        ValueError: While trying to set the hidden weights of the artificial \
neural network `ann` of element `?`, the following error occured: could not \
broadcast input array from shape (3,2) into shape (2,3,3)
        """
        return numpy.asarray(self._cann.weights_hidden)

    def _set_weights_hidden(self, values):
        if values is None:
            self._del_weights_hidden()
        else:
            try:
                self._cann.weights_hidden = numpy.full(
                    self.shape_weights_hidden, values, dtype=float)
            except BaseException:
                objecttools.augment_excmessage(
                    'While trying to set the hidden weights of the artificial '
                    'neural network `%s` of element `%s`'
                    % (self.name, objecttools.devicename(self)))

    def _del_weights_hidden(self):
        self._cann.weights_hidden = numpy.full(self.shape_weights_hidden,
                                               numpy.nan)
        for idx_layer in range(self.nmb_layers-1):
            for idx_neuron1 in range(self.nmb_neurons[idx_layer]):
                for idx_neuron2 in range(self.nmb_neurons[idx_layer+1]):
                    self._cann.weights_hidden[idx_layer,
                                              idx_neuron1,
                                              idx_neuron2] = 0.

    weights_hidden = exceptiontools.dependent_property(
        'weights_hidden', _get_weights_hidden,
        _set_weights_hidden, _del_weights_hidden)

    @property
    def shape_weights_hidden(self):
        """Shape of the array containing the activation of the hidden neurons.

        >>> from hydpy import dummies
        >>> dummies.ann.shape_weights_hidden
        (1, 2, 2)
        """
        if self.nmb_layers > 1:
            return (self.nmb_layers-1,
                    self._max_nmb_neurons,
                    self._max_nmb_neurons)
        return (0, 0, 0)

    @property
    def nmb_weights_hidden(self):
        """Number of hidden weights.

        >>> from hydpy import dummies
        >>> dummies.ann.nmb_weights_hidden
        2
        """
        nmb = 0
        for idx_layer in range(self.nmb_layers-1):
            nmb += self.nmb_neurons[idx_layer] * self.nmb_neurons[idx_layer+1]
        return nmb

    def _get_intercepts_hidden(self):
        """Intercepts of all neurons of the hidden layers.

        All intercepts are handled in a 1-dimensional array:

        >>> from hydpy import ANN
        >>> ann = ANN()
        >>> ann(nmb_neurons=(3, 2))
        >>> ann.intercepts_hidden
        array([[  0.,   0.,   0.],
               [  0.,   0.,  nan]])

        It is allowed to set values via slicing:

        >>> ann.intercepts_hidden[0, :] = 1.
        >>> ann.intercepts_hidden
        array([[  1.,   1.,   1.],
               [  0.,   0.,  nan]])

        If possible, type conversions are performed:

        >>> ann.intercepts_hidden = '2'
        >>> ann.intercepts_hidden
        array([[ 2.,  2.,  2.],
               [ 2.,  2.,  2.]])

        One can assign whole matrices directly:

        >>> import numpy
        >>> ann.intercepts_hidden = [1.0, 3.0, 2.0]
        >>> ann.intercepts_hidden
        array([[ 1.,  3.,  2.],
               [ 1.,  3.,  2.]])

        One can also delete the values contained in the array:

        >>> del ann.intercepts_hidden
        >>> ann.intercepts_hidden
        array([[  0.,   0.,   0.],
               [  0.,   0.,  nan]])

        Errors like wrong shapes (or unconvertible inputs) result in error
        messages:

        >>> ann.intercepts_hidden = [1.0, 3.0]
        Traceback (most recent call last):
        ...
        ValueError: While trying to set the neuron related intercepts of the \
artificial neural network `ann` of element `?`, the following error occured: \
could not broadcast input array from shape (2) into shape (2,3)

        The number of input intercepts is available as a property:

        >>> ann.nmb_intercepts_hidden
        5
        """
        return numpy.asarray(self._cann.intercepts_hidden)

    def _set_intercepts_hidden(self, values):
        if values is None:
            self._del_intercepts_hidden()
        else:
            try:
                self._cann.intercepts_hidden = numpy.full(
                    self.shape_intercepts_hidden, values, dtype=float)
            except BaseException:
                objecttools.augment_excmessage(
                    'While trying to set the neuron related intercepts of '
                    'the artificial neural network `%s` of element `%s`'
                    % (self.name, objecttools.devicename(self)))

    def _del_intercepts_hidden(self):
        self._cann.intercepts_hidden = numpy.full(
            self.shape_intercepts_hidden,
            numpy.nan)
        for idx_layer in range(self.nmb_layers):
            for idx_neuron in range(self.nmb_neurons[idx_layer]):
                self._cann.intercepts_hidden[idx_layer, idx_neuron] = 0.

    intercepts_hidden = exceptiontools.dependent_property(
        'intercepts_hidden', _get_intercepts_hidden,
        _set_intercepts_hidden, _del_intercepts_hidden)

    @property
    def shape_intercepts_hidden(self):
        """Shape if the array containing the intercepts of neurons of
        the hidden layers."""
        return (self.nmb_layers, self._max_nmb_neurons)

    @property
    def nmb_intercepts_hidden(self):
        """Number of input intercepts."""
        return sum(self.nmb_neurons)

    def _get_intercepts_output(self):
        """Intercepts of all output nodes.

        All intercepts are handled in a 1-dimensional array:

        >>> from hydpy import ANN
        >>> ann = ANN()
        >>> ann(nmb_outputs=3)
        >>> ann.intercepts_output
        array([ 0.,  0.,  0.])

        It is allowed to set values via slicing:

        >>> ann.intercepts_output[1:] = 1.
        >>> ann.intercepts_output
        array([ 0.,  1.,  1.])

        If possible, type conversions are performed:

        >>> ann.intercepts_output = '2'
        >>> ann.intercepts_output
        array([ 2.,  2.,  2.])

        One can assign whole matrices directly:

        >>> import numpy
        >>> ann.intercepts_output = [1.0, 3.0, 2.0]
        >>> ann.intercepts_output
        array([ 1.,  3.,  2.])

        One can also delete the values contained in the array:

        >>> del ann.intercepts_output
        >>> ann.intercepts_output
        array([ 0.,  0.,  0.])

        Errors like wrong shapes (or unconvertible inputs) result in error
        messages:

        >>> ann.intercepts_output = [1.0, 3.0]
        Traceback (most recent call last):
        ...
        ValueError: While trying to set the output node related intercepts \
of the artificial neural network `ann` of element `?`, the following error \
occured: could not broadcast input array from shape (2) into shape (3)
        """
        return numpy.asarray(self._cann.intercepts_output)

    def _set_intercepts_output(self, values):
        if values is None:
            self._del_intercepts_output()
        else:
            try:
                self._cann.intercepts_output = numpy.full(
                    self.shape_intercepts_output, values, dtype=float)
            except BaseException:
                objecttools.augment_excmessage(
                    'While trying to set the output node related intercepts '
                    'of the artificial neural network `%s` of element `%s`'
                    % (self.name, objecttools.devicename(self)))

    def _del_intercepts_output(self):
        self._cann.intercepts_output = numpy.zeros(
            self.shape_intercepts_output)

    intercepts_output = exceptiontools.dependent_property(
        'intercepts_output', _get_intercepts_output,
        _set_intercepts_output, _del_intercepts_output)

    @property
    def shape_intercepts_output(self):
        """Shape if the array containing the intercepts of neurons of
        the hidden layers.

        >>> from hydpy import dummies
        >>> dummies.ann.shape_intercepts_output
        (1,)
        """
        return (self.nmb_outputs,)

    @property
    def nmb_intercepts_output(self):
        """Number of output intercepts.

        >>> from hydpy import dummies
        >>> dummies.ann.nmb_intercepts_output
        1
        """
        return self.nmb_outputs

    def _get_inputs(self):
        """Values of the input nodes.

        All input values are handled in a 1-dimensional array:

        >>> from hydpy import ANN
        >>> ann = ANN()
        >>> ann(nmb_inputs=3)
        >>> ann.inputs
        array([ 0.,  0.,  0.])

        It is allowed to set values via slicing:

        >>> ann.inputs[1:] = 1.
        >>> ann.inputs
        array([ 0.,  1.,  1.])

        If possible, type conversions are performed:

        >>> ann.inputs = '2'
        >>> ann.inputs
        array([ 2.,  2.,  2.])

        One can assign whole matrices directly:

        >>> import numpy
        >>> ann.inputs = [1.0, 3.0, 2.0]
        >>> ann.inputs
        array([ 1.,  3.,  2.])

        One can also delete the values contained in the array:

        >>> del ann.inputs
        >>> ann.inputs
        array([ 0.,  0.,  0.])

        Errors like wrong shapes (or unconvertible inputs) result in error
        messages:

        >>> ann.inputs = [1.0, 3.0]
        Traceback (most recent call last):
        ...
        ValueError: While trying to set the inputs of the artificial neural \
network `ann` of element `?`, the following error occured: could not \
broadcast input array from shape (2) into shape (3)
        """
        return numpy.asarray(self._cann.inputs)

    def _set_inputs(self, values):
        try:
            self._cann.inputs = numpy.full(self.nmb_inputs,
                                           values, dtype=float)
        except BaseException:
            objecttools.augment_excmessage(
                'While trying to set the inputs of the artificial '
                'neural network `%s` of element `%s`'
                % (self.name, objecttools.devicename(self)))

    def _del_inputs(self):
        self._cann.inputs = numpy.zeros(self.nmb_inputs)

    inputs = exceptiontools.dependent_property(
        'inputs', _get_inputs, _set_inputs, _del_inputs)

    def _get_outputs(self):
        """Values of the output nodes.

        All output values are handled in a 1-dimensional array:

        >>> from hydpy import ANN
        >>> ann = ANN()
        >>> ann(nmb_outputs=3)
        >>> ann.outputs
        array([ 0.,  0.,  0.])

        It is not allowed to change output values manually:

        >>> ann.outputs = 1.0
        Traceback (most recent call last):
        ...
        AttributeError: Attribute `outputs` of object `ann` \
cannot be used this way.
        """
        return numpy.asarray(self._cann.outputs)

    def _del_outputs(self):
        self._cann.outputs = numpy.zeros(self.nmb_outputs)

    outputs = exceptiontools.dependent_property(
        'outputs', _get_outputs, fdel=_del_outputs)

    def _get_neurons(self):
        """The activation of the neurons of the hidden layers.

        >>> from hydpy import dummies
        >>> dummies.ann.neurons
        array([[ 1.,  1.],
               [ 0.,  0.]])
        """
        return numpy.array(self._cann.neurons)

    def _del_neurons(self):
        nmb_neurons = numpy.asarray(self._cann.nmb_neurons)
        self._cann.neurons = numpy.zeros((self.nmb_layers, max(nmb_neurons)))

    neurons = exceptiontools.dependent_property(
        'neurons', _get_neurons, fdel=_del_neurons)

    def process_actual_input(self):
        """Calculates the network output values based on the input values
        defined previously.

        For more information see the documentation on class |anntools.ANN|.
        """
        self._cann.process_actual_input()

    @property
    def nmb_weights(self):
        """Number of all input, inner, and output weights.

        >>> from hydpy import dummies
        >>> dummies.ann.nmb_weights
        7
        """
        return (self.nmb_weights_input +
                self.nmb_weights_hidden +
                self.nmb_weights_output)

    @property
    def nmb_intercepts(self):
        """Number of all inner and output intercepts.

        >>> from hydpy import dummies
        >>> dummies.ann.nmb_intercepts
        4
        """
        return self.nmb_intercepts_hidden + self.nmb_intercepts_output

    @property
    def nmb_parameters(self):
        """Sum of |anntools.ANN.nmb_weights| and |anntools.ANN.nmb_intercepts|.

        >>> from hydpy import dummies
        >>> dummies.ann.nmb_parameters
        11
        """
        return self.nmb_weights + self.nmb_intercepts

    def verify(self):
        """Raise a |RuntimeError| if the network's shape is not defined
        completely.


        >>> from hydpy import dummies
        >>> dummies.ann.verify()

        >>> from hydpy import ANN
        >>> ANN().verify()
        Traceback (most recent call last):
        ...
        RuntimeError: The shape of the the artificial neural network \
parameter `ann` of element `?` has not been defined so far.
        """
        if not self._isready:
            raise RuntimeError(
                'The shape of the the artificial neural network '
                'parameter `%s`%shas not been defined so far.'
                % (self.name, objecttools.elementphrase(self)))

    def assignrepr(self, prefix):
        """Return a string representation of the actual |anntools.ANN| object
        that is prefixed with the given string."""
        prefix = '%s%s(' % (prefix, self.name)
        blanks = len(prefix)*' '
        lines = [(objecttools.assignrepr_value(
            self.nmb_inputs, '%snmb_inputs=' % prefix)+',')]
        lines.append(objecttools.assignrepr_tuple(
            self.nmb_neurons, '%snmb_neurons=' % blanks)+',')
        lines.append(objecttools.assignrepr_value(
            self.nmb_outputs, '%snmb_outputs=' % blanks)+',')
        lines.append(objecttools.assignrepr_list2(
            self.weights_input, '%sweights_input=' % blanks)+',')
        if self.nmb_layers > 1:
            lines.append(objecttools.assignrepr_list3(
                self.weights_hidden, '%sweights_hidden=' % blanks)+',')
        lines.append(objecttools.assignrepr_list2(
            self.weights_output, '%sweights_output=' % blanks)+',')
        lines.append(objecttools.assignrepr_list2(
            self.intercepts_hidden, '%sintercepts_hidden=' % blanks)+',')
        lines.append(objecttools.assignrepr_list(
            self.intercepts_output, '%sintercepts_output=' % blanks)+')')
        return '\n'.join(lines)

    def __repr__(self):
        return self.assignrepr(prefix='')

    def plot(self, xmin, xmax, idx_input=0, idx_output=0, points=100,
             **kwargs):
        """Plot the relationship between a certain input (`idx_input`) and a
        certain output (`idx_output`) variable described by the actual
        |anntools.ANN| object.

        Define the lower and the upper bound of the x axis via arguments
        `xmin` and `xmax`.  The number of plotting points can be modified
        by argument `points`.  Additional `matplotlib` plotting arguments
        can be passed as keyword arguments.
        """
        xs_ = numpy.linspace(xmin, xmax, points)
        ys_ = numpy.zeros(xs_.shape)
        for idx, x__ in enumerate(xs_):
            self.inputs[idx_input] = x__
            self.process_actual_input()
            ys_[idx] = self.outputs[idx_output]
        pyplot.plot(xs_, ys_, **kwargs)


abctools.ParameterABC.register(ANN)
abctools.ANNABC.register(ANN)


def ann(**kwargs):
    """Return a new stand alone |anntools.ANN| object with the given parameter
    values.

    The purpose of this function is to allow for string representations of
    parameters containing multiple |anntools.ANN| instances.

    When passing no arguments, the default values of class |anntools.ANN| will
    be applied:

    >>> from hydpy import ANN
    >>> ann1 = ann()
    >>> ann1
    ann(nmb_inputs=1,
        nmb_neurons=(1,),
        nmb_outputs=1,
        weights_input=[[0.0]],
        weights_output=[[0.0]],
        intercepts_hidden=[[0.0]],
        intercepts_output=[0.0])

    Of course, all parameter values can be changed:

    >>> ann2 = ann(nmb_inputs=1, nmb_neurons=(1,), nmb_outputs=1,
    ...            weights_input=4.0, weights_output=3.0,
    ...            intercepts_hidden=-16.0, intercepts_output=-1.0)
    >>> ann2
    ann(nmb_inputs=1,
        nmb_neurons=(1,),
        nmb_outputs=1,
        weights_input=[[4.0]],
        weights_output=[[3.0]],
        intercepts_hidden=[[-16.0]],
        intercepts_output=[-1.0])

    The following line is just thought to make clear, that two independent
    |anntools.ANN| objects have been initialized (instead of changing the
    values of an existing |anntools.ANN| object vai its `call` method):

    >>> ann1 is ann2
    False
    """
    new_ann = ANN()
    new_ann(**kwargs)
    return new_ann


class SeasonalANN(object):
    """Handles relationships described by artificial neural networks that
    vary within an anual cycle.

    Class |anntools.SeasonalANN| is an alternative implementation of class
    |SeasonalParameter| specifically designed for handling multiple
    |anntools.ANN| objects that are valid for different times of the year,
    described by |TOY| objects.  The total output of a |anntools.SeasonalANN|
    object is a weighted mean of the output of one or two "normal" neural
    networks.  |anntools.SeasonalANN.ratios| used for weighting depend
    on the actual time of the year.

    To explain this in more detail, let us define a |anntools.SeasonalANN|
    object first, that contains three "normal" networks for January, 1,
    March, 1, and July, 1, respectively (note that this example is similar
    to the example used to describe class
    |SeasonalParameter|):

    >>> from hydpy import SeasonalANN, ann
    >>> seasonalann = SeasonalANN()
    >>> seasonalann.simulationstep = '1d'
    >>> seasonalann(
    ...     _1_1_12=ann(nmb_inputs=1, nmb_neurons=(1,), nmb_outputs=1,
    ...                 weights_input=0.0, weights_output=0.0,
    ...                 intercepts_hidden=0.0, intercepts_output=1.0),
    ...     _7_1_12=ann(nmb_inputs=1, nmb_neurons=(1,), nmb_outputs=1,
    ...                 weights_input=4.0, weights_output=3.0,
    ...                 intercepts_hidden=-16.0, intercepts_output=-1.0),
    ...     _3_1_12=ann(nmb_inputs=1, nmb_neurons=(1,), nmb_outputs=1,
    ...                 weights_input=0.0, weights_output=0.0,
    ...                 intercepts_hidden=0.0, intercepts_output=-1.0))

    The confused time order in the initialization call above does not pose
    a problem, as |anntools.SeasonalANN| performs time sorting internally:

    >>> seasonalann
    seasonalann(toy_1_1_12_0_0=ann(nmb_inputs=1,
                                   nmb_neurons=(1,),
                                   nmb_outputs=1,
                                   weights_input=[[0.0]],
                                   weights_output=[[0.0]],
                                   intercepts_hidden=[[0.0]],
                                   intercepts_output=[1.0]),
                toy_3_1_12_0_0=ann(nmb_inputs=1,
                                   nmb_neurons=(1,),
                                   nmb_outputs=1,
                                   weights_input=[[0.0]],
                                   weights_output=[[0.0]],
                                   intercepts_hidden=[[0.0]],
                                   intercepts_output=[-1.0]),
                toy_7_1_12_0_0=ann(nmb_inputs=1,
                                   nmb_neurons=(1,),
                                   nmb_outputs=1,
                                   weights_input=[[4.0]],
                                   weights_output=[[3.0]],
                                   intercepts_hidden=[[-16.0]],
                                   intercepts_output=[-1.0]))

    The property |anntools.SeasonalANN.shape| does reflect the number of
    required weighting ratios for each time of year (in this example:
    366 days per year) and each neural network (in this example: three):

    >>> seasonalann.shape
    (366, 3)

    For safety reasons, |anntools.SeasonalANN.shape| should normally not
    be changed manually:

    >>> seasonalann.shape = (366, 4)
    Traceback (most recent call last):
    ...
    AttributeError: can't set attribute

    The following interactive shows how the |anntools.SeasonalANN.ratios|
    used for weighting are calculated:

    .. testsetup::

        >>> from bokeh import plotting, models, palettes
        ...
        >>> from hydpy import docs
        >>> import os
        >>> plotting.output_file(os.path.join(
        ...     docs.__path__[0], 'html', 'anntools.SeasonalANN.ratios.html'))
        >>> hover = models.HoverTool(tooltips=[
        ...     ("(x,y)", "($x, $y)")])
        >>> plot = plotting.figure(toolbar_location="above",
        ...                        plot_width=500, plot_height=300)
        >>> plot.tools.append(hover)
        >>> legend_entries = []
        >>> for idx, (toy, color) in enumerate(
        ...         zip(seasonalann.toys, palettes.Dark2_5)):
        ...     line = plot.line(range(366), seasonalann.ratios[:, idx],
        ...                      alpha=0.8, muted_alpha=0.2, color=color)
        ...     line.muted = True
        ...     legend_entries.append((str(toy), [line]))
        >>> legend = models.Legend(items=legend_entries,
        ...                        location=(10, 0),
        ...                        click_policy='mute')
        >>> plot.add_layout(legend, 'right')
        >>> label_dict = {0: 'Jan 1',
        ...               60: 'Mar 1',
        ...               182: 'Jul 1'}
        >>> plot.xaxis.ticker =  sorted(label_dict.keys())
        >>> plot.xaxis.formatter = models.FuncTickFormatter(
        ...     code='var labels = %s; return labels[tick];' % label_dict)
        >>> dummy = plotting.save(plot)

    .. raw:: html

        <iframe
            src="anntools.SeasonalANN.ratios.html"
            width="100%"
            height="300px"
            frameborder=0
        ></iframe>

    For example, on July, 1 (which is the 183th day of a leap year),
    only the output of the third network is relevant:

    >>> from hydpy import print_values
    >>> print_values(seasonalann.ratios[182])
    0.0, 0.0, 1.0

    On Juni, 30, and July, 2, also the second and the first neural network
    are relevant, respectively:

    >>> print_values(seasonalann.ratios[181])
    0.0, 0.008197, 0.991803
    >>> print_values(seasonalann.ratios[183])
    0.005435, 0.0, 0.994565

    Inserting data, processing this data, and fetching the output works
    as explained for class |anntools.ANN|, except that the index of the
    actual time of year needs to be passed as the single argument of
    |anntools.SeasonalANN.process_actual_input|.  Passing the index value
    `182` activates the third network only, which is configured exactly
    as the one exemplifying class |anntools.ANN|:

    >>> from hydpy import round_
    >>> for input_ in range(9):
    ...     seasonalann.inputs[0] = input_
    ...     seasonalann.process_actual_input(182)
    ...     round_([input_, seasonalann.outputs[0]])
    0, -1.0
    1, -0.999982
    2, -0.998994
    3, -0.946041
    4, 0.5
    5, 1.946041
    6, 1.998994
    7, 1.999982
    8, 2.0

    To see that the final output values are actually the weighted mean
    of the output values of the single neural networks, we repeat the
    above example for January, 13, where the first and the second neural
    network have ratios of 0.8 and 0.2 respectively:

    >>> print_values(seasonalann.ratios[12])
    0.8, 0.2, 0.0

    For both networks all parameters except the output intercepts are
    zero.  Hence, the calculated output is independent of the given input.
    The output of the first network (1.0) dominates the output of the
    second network (-1.0):

    >>> from hydpy import round_
    >>> for input_ in range(9):
    ...     seasonalann.inputs[0] = input_
    ...     seasonalann.process_actual_input(12)
    ...     round_([input_, seasonalann.outputs[0]])
    0, 0.6
    1, 0.6
    2, 0.6
    3, 0.6
    4, 0.6
    5, 0.6
    6, 0.6
    7, 0.6
    8, 0.6

    It is of great importance that all contained neural networks are
    consistent.  Hence some tests are performed:

    >>> seasonalann = SeasonalANN()
    >>> seasonalann.process_actual_input(0)
    Traceback (most recent call last):
    ...
    RuntimeError: The seasonal neural network collection `seasonalann` \
of element `?` has not been properly prepared so far.

    >>> seasonalann(1)
    Traceback (most recent call last):
    ...
    TypeError: Type `int` is not (a subclass of) type `ANN`.

    >>> seasonalann(
    ...     _13_1_12=ann(nmb_inputs=2, nmb_neurons=(1,), nmb_outputs=1,
    ...                  weights_input=0.0, weights_output=0.0,
    ...                  intercepts_hidden=0.0, intercepts_output=1.0))
    Traceback (most recent call last):
    ...
    ValueError: While trying to add a season specific neural network to \
parameter `seasonalann` of element `?`, the following error occured: \
While trying to retrieve the month for TOY (time of year) object based \
on the string `_13_1_12`, the following error occured: \
The value of property `month` of TOY (time of year) objects must lie \
within the range `(1, 12)`, but the given value is `13`.

    >>> seasonalann(
    ...     ann(nmb_inputs=2, nmb_neurons=(1,), nmb_outputs=1,
    ...         weights_input=0.0, weights_output=0.0,
    ...         intercepts_hidden=0.0, intercepts_output=1.0))
    >>> seasonalann
    seasonalann(ann(nmb_inputs=2,
                    nmb_neurons=(1,),
                    nmb_outputs=1,
                    weights_input=[[0.0],
                                   [0.0]],
                    weights_output=[[0.0]],
                    intercepts_hidden=[[0.0]],
                    intercepts_output=[1.0]))

    >>> seasonalann(
    ...     ann(nmb_inputs=2, nmb_neurons=(1,), nmb_outputs=1,
    ...         weights_input=0.0, weights_output=0.0,
    ...         intercepts_hidden=0.0, intercepts_output=1.0),
    ...     _7_1_12=ann(nmb_inputs=1, nmb_neurons=(1,), nmb_outputs=1,
    ...                 weights_input=4.0, weights_output=3.0,
    ...                 intercepts_hidden=-16.0, intercepts_output=-1.0),
    ...     _3_1_12=ann(nmb_inputs=1, nmb_neurons=(1,), nmb_outputs=1,
    ...                 weights_input=0.0, weights_output=0.0,
    ...                 intercepts_hidden=0.0, intercepts_output=-1.0))
    Traceback (most recent call last):
    ...
    ValueError: Type `SeasonalANN` accepts either a single positional \
argument or an arbitrary number of keyword arguments, but for the \
corresponding parameter of element `?` 1 positional and 2 keyword \
arguments have been given.

    >>> seasonalann(
    ...     _1_1_12=ann(nmb_inputs=2, nmb_neurons=(1,), nmb_outputs=1,
    ...                 weights_input=0.0, weights_output=0.0,
    ...                 intercepts_hidden=0.0, intercepts_output=1.0),
    ...     _7_1_12=ann(nmb_inputs=1, nmb_neurons=(1,), nmb_outputs=1,
    ...                 weights_input=4.0, weights_output=3.0,
    ...                 intercepts_hidden=-16.0, intercepts_output=-1.0),
    ...     _3_1_12=ann(nmb_inputs=1, nmb_neurons=(1,), nmb_outputs=1,
    ...                 weights_input=0.0, weights_output=0.0,
    ...                 intercepts_hidden=0.0, intercepts_output=-1.0))
    Traceback (most recent call last):
    ...
    RuntimeError: The number of input and output values of all neural \
networks contained by a seasonal neural network collection must be \
identical and be known by the containing object.  But the seasonal \
neural network collection `seasonalann` of element `?` assumes `2` input \
and `1` output values, while the network corresponding to the time of \
year `toy_3_1_12_0_0` requires `1` input and `1` output values.

    Whenever a test fails, all networks are removed for safety:

    >>> seasonalann
    seasonalann()

    Alternatively, neural networks can be added individually via
    attribute access:

    >>> jan = ann(nmb_inputs=1, nmb_neurons=(1,), nmb_outputs=1,
    ...           weights_input=0.0, weights_output=0.0,
    ...           intercepts_hidden=0.0, intercepts_output=1.0)
    >>> seasonalann.toy_1_1_12 = jan

    Setting an attribute updates everything, e.g.:

    >>> round_(seasonalann.ratios[0])
    1.0

    The mentioned safety checks do also apply when adding networks
    via attribute access, e.g.:

    >>> seasonalann.toy_7_1_12 = ann(nmb_inputs=2,
    ...                              nmb_neurons=(1,),
    ...                              nmb_outputs=1,
    ...                              weights_input=0.0,
    ...                              weights_output=0.0,
    ...                              intercepts_hidden=0.0,
    ...                              intercepts_output=1.0)
    Traceback (most recent call last):
    ...
    RuntimeError: While trying to assign a new neural network to the \
seasonal neural network collection `seasonalann` of element `?` based \
on name `toy_7_1_12`, the following error occured: \
The number of input and output values of all neural networks contained \
by a seasonal neural network collection must be identical and be known \
by the containing object.  But the seasonal neural network collection \
`seasonalann` of element `?` assumes `1` input and `1` output values, \
while the network corresponding to the time of year `toy_7_1_12_0_0` \
requires `2` input and `1` output values.

    Besides setting new networks, getting and deleting them are also
    suppported:

    >>> seasonalann.toy_1_1_12 = jan
    >>> seasonalann.toy_1_1_12
    ann(nmb_inputs=1,
        nmb_neurons=(1,),
        nmb_outputs=1,
        weights_input=[[0.0]],
        weights_output=[[0.0]],
        intercepts_hidden=[[0.0]],
        intercepts_output=[1.0])
    >>> del seasonalann.toy_1_1_12

    These error messages related to attribute access are provided:

    >>> seasonalann.toy_1_1_12
    Traceback (most recent call last):
    ...
    AttributeError: While trying to look up for a neural network handled \
by the seasonal neural network collection `seasonalann` of element `?` \
based on name `toy_1_1_12`, the following error occured: No neural network \
is registered under a TOY object named `toy_1_1_12_0_0`.

    >>> del seasonalann.toy_1_1_12
    Traceback (most recent call last):
    ...
    AttributeError: While trying to remove a new neural network from the \
seasonal neural network collection `seasonalann` of element `?` based on \
name `toy_1_1_12`, the following error occured: No neural network is \
registered under a TOY object named `toy_1_1_12_0_0`.

    >>> seasonalann.toy_1_1_12 = 1
    Traceback (most recent call last):
    ...
    TypeError: While trying to assign a new neural network to the seasonal \
neural network collection `seasonalann` of element `?` based on name \
`toy_1_1_12`, the following error occured: Value `1` of type `int` has \
been given, but a value of type `ANN` is required.

    Setting and deleting "normal" attributes is supported:

    >>> seasonalann.temp = 999
    >>> seasonalann.temp
    999
    >>> del seasonalann.temp
    >>> seasonalann.temp
    Traceback (most recent call last):
    ...
    AttributeError: 'SeasonalANN' object has no attribute 'temp'
    """
    NDIM = 0
    TYPE = 'annutils.SeasonalANN'
    TIME = None
    SPAN = (None, None)

    parameterstep = parametertools.Parameter.__dict__['parameterstep']
    simulationstep = parametertools.Parameter.__dict__['simulationstep']

    def __init__(self):
        self.subpars = None
        self.fastaccess = objecttools.FastAccess()
        self._toy2ann = {}
        self.__sann = None
        self._do_refresh = True

    def connect(self, subpars):
        """Connect the actual |anntools.SeasonalANN| object with the given
        |SubParameters| object."""
        self.subpars = subpars
        self.fastaccess = subpars.fastaccess

    name = property(objecttools.name)

    def __call__(self, *args, **kwargs):
        self._toy2ann.clear()
        self._do_refresh = False
        try:
            if (len(args) > 1) or (args and kwargs):
                raise ValueError(
                    'Type `%s` accepts either a single positional argument or '
                    'an arbitrary number of keyword arguments, but for the '
                    'corresponding parameter of element `%s` %d positional '
                    'and %d keyword arguments have been given.'
                    % (objecttools.classname(self),
                       objecttools.devicename(self),
                       len(args), len(kwargs)))
            if args:
                kwargs['_1'] = args[0]
            for (toystr, value) in kwargs.items():
                if not isinstance(value, abctools.ANNABC):
                    raise TypeError(
                        'Type `%s` is not (a subclass of) type `ANN`.'
                        % objecttools.classname(value))
                try:
                    setattr(self, str(timetools.TOY(toystr)), value)
                except BaseException:
                    objecttools.augment_excmessage(
                        'While trying to add a season specific neural '
                        'network to parameter `%s` of element `%s`'
                        % (self.name, objecttools.devicename(self)))
        except BaseException as exc:
            self._toy2ann.clear()
            raise exc
        finally:
            self._do_refresh = True
            self.refresh()

    def refresh(self):
        """Prepare the actual |anntools.SeasonalANN| object for calculations.

        Dispite all automated refreshings explained in the general
        documentation on class |anntools.SeasonalANN|, it is still possible
        to destroy the inner consistency of a |anntools.SeasonalANN| instance,
        as it stores its |anntools.ANN| objects by reference.  This is shown
        by the following example:

        >>> from hydpy import SeasonalANN, ann
        >>> seasonalann = SeasonalANN()
        >>> seasonalann.simulationstep = '1d'
        >>> jan = ann(nmb_inputs=1, nmb_neurons=(1,), nmb_outputs=1,
        ...           weights_input=0.0, weights_output=0.0,
        ...           intercepts_hidden=0.0, intercepts_output=1.0)
        >>> seasonalann(_1_1_12=jan)
        >>> jan.nmb_inputs, jan.nmb_outputs = 2, 3
        >>> jan.nmb_inputs, jan.nmb_outputs
        (2, 3)
        >>> seasonalann.nmb_inputs, seasonalann.nmb_outputs
        (1, 1)

        Due to the C level implementation of the mathematical core of
        both |anntools.ANN| and |anntools.SeasonalANN| in module |annutils|,
        such an inconsistency might result in a program crash without any
        informative error message.  Whenever you are afraid some
        inconsistency might have crept in, and you want to repair it,
        call method |anntools.SeasonalANN.refresh| explicitly:

        >>> seasonalann.refresh()
        >>> jan.nmb_inputs, jan.nmb_outputs
        (2, 3)
        >>> seasonalann.nmb_inputs, seasonalann.nmb_outputs
        (2, 3)
        """
        if self._do_refresh:
            if self.anns:
                self.__sann = annutils.SeasonalANN(self.anns)
                setattr(self.fastaccess, self.name, self._sann)
                self._setshape((None, self._sann.nmb_anns))
                if self._sann.nmb_anns > 1:
                    self._interp()
                else:
                    self._sann.ratios[:, 0] = 1.
                self.verify()
            else:
                self.__sann = None

    def verify(self):
        """Raise a |RuntimeError| and removes all handled neural networks,
        if the they are defined inconsistently.

        Dispite all automated safety checks explained in the general
        documentation on class |anntools.SeasonalANN|, it is still possible
        to destroy the inner consistency of a |anntools.SeasonalANN| instance,
        as it stores its |anntools.ANN| objects by reference.  This is shown
        by the following example:

        >>> from hydpy import SeasonalANN, ann
        >>> seasonalann = SeasonalANN()
        >>> seasonalann.simulationstep = '1d'
        >>> jan = ann(nmb_inputs=1, nmb_neurons=(1,), nmb_outputs=1,
        ...           weights_input=0.0, weights_output=0.0,
        ...           intercepts_hidden=0.0, intercepts_output=1.0)
        >>> seasonalann(_1_1_12=jan)
        >>> jan.nmb_inputs, jan.nmb_outputs = 2, 3
        >>> jan.nmb_inputs, jan.nmb_outputs
        (2, 3)
        >>> seasonalann.nmb_inputs, seasonalann.nmb_outputs
        (1, 1)

        Due to the C level implementation of the mathematical core of both
        |anntools.ANN| and |anntools.SeasonalANN| in module |annutils|,
        such an inconsistency might result in a program crash without any
        informative error message.  Whenever you are afraid some
        inconsistency might have crept in, and you want to find out if this
        is actually the case, call method |anntools.SeasonalANN.verify|
        explicitly:

        >>> seasonalann.verify()
        Traceback (most recent call last):
        ...
        RuntimeError: The number of input and output values of all neural \
networks contained by a seasonal neural network collection must be \
identical and be known by the containing object.  But the seasonal \
neural network collection `seasonalann` of element `?` assumes `1` input \
and `1` output values, while the network corresponding to the time of \
year `toy_1_1_12_0_0` requires `2` input and `3` output values.

        >>> seasonalann
        seasonalann()

        >>> seasonalann.verify()
        Traceback (most recent call last):
        ...
        RuntimeError: Seasonal artificial neural network collections need \
to handle at least one "normal" single neural network, but for the seasonal \
neural network `seasonalann` of element `?` none has been defined so far.

        """
        if not self.anns:
            self._toy2ann.clear()
            raise RuntimeError(
                'Seasonal artificial neural network collections need '
                'to handle at least one "normal" single neural network, '
                'but for the seasonal neural network `%s` of element '
                '`%s` none has been defined so far.'
                % (self.name, objecttools.devicename(self)))
        for toy, ann_ in self:
            ann_.verify()
            if ((self.nmb_inputs != ann_.nmb_inputs) or
                    (self.nmb_outputs != ann_.nmb_outputs)):
                self._toy2ann.clear()
                raise RuntimeError(
                    'The number of input and output values of all neural '
                    'networks contained by a seasonal neural network '
                    'collection must be identical and be known by the '
                    'containing object.  But the seasonal neural '
                    'network collection `%s` of element `%s` assumes '
                    '`%d` input and `%d` output values, while the network '
                    'corresponding to the time of year `%s` requires '
                    '`%d` input and `%d` output values.'
                    % (self.name, objecttools.devicename(self),
                       self.nmb_inputs, self.nmb_outputs,
                       toy,
                       ann_.nmb_inputs, ann_.nmb_outputs))

    def _interp(self):
        ratios = self.ratios
        ratios[:, :] = 0.0
        toys = self.toys
        timegrid = timetools.Timegrid(
            timetools.TOY._STARTDATE+self.simulationstep/2,
            timetools.TOY._ENDDATE+self.simulationstep/2,
            self.simulationstep)
        for tdx, date in enumerate(timegrid):
            xnew = timetools.TOY(date)
            for idx_1, x_1 in enumerate(toys):
                if x_1 > xnew:
                    idx_0 = idx_1-1
                    x_0 = toys[idx_0]
                    break
            else:
                idx_0 = -1
                idx_1 = 0
                x_0 = toys[idx_0]
                x_1 = toys[idx_1]
            ratios[tdx, idx_1] = (xnew-x_0)/(x_1-x_0)
            ratios[tdx, idx_0] = 1.-ratios[tdx, idx_1]

    def _getshape(self):
        return tuple(int(sub) for sub in self.ratios.shape)

    def _setshape(self, shape):
        try:
            shape = (int(shape),)
        except TypeError:
            pass
        shp = list(shape)
        shp[0] = timetools.Period('366d')/self.simulationstep
        shp[0] = int(numpy.ceil(round(shp[0], 10)))
        getattr(self.fastaccess, self.name).ratios = numpy.zeros(
            shp, dtype=float)

    shape = property(
        _getshape,
        doc='The shape of array |anntools.SeasonalANN.ratios|.')

    @property
    def toys(self):
        """A sorted |tuple| of all contained |TOY| objects."""
        return tuple(toy for (toy, ann) in self)

    @property
    def anns(self):
        """A sorted |tuple| of all contained |anntools.ANN| objects."""
        return tuple(ann for (toy, ann) in self)

    @property
    def ratios(self):
        """Ratios for weighting the single neural network outputs."""
        return numpy.asarray(self._sann.ratios)

    @property
    def _sann(self):
        sann = self.__sann
        if sann:
            return sann
        else:
            raise RuntimeError(
                'The seasonal neural network collection `%s` of '
                'element `%s` has not been properly prepared so far.'
                % (self.name, objecttools.devicename(self)))

    @property
    def nmb_inputs(self):
        """Number of input values of all neural networks."""
        return self._sann.nmb_inputs

    @property
    def inputs(self):
        """General input data for all neural networks."""
        return numpy.asarray(self._sann.inputs)

    @property
    def nmb_outputs(self):
        """Number of output values of all neural networks."""
        return self._sann.nmb_outputs

    @property
    def outputs(self):
        """Weighted output of the individual neural networks."""
        return numpy.asarray(self._sann.outputs)

    def process_actual_input(self, idx_toy):
        """Calculate the network output values based on the input values
        defined previously for the given index referencing the actual
        time of year.
        """
        self._sann.process_actual_input(idx_toy)

    def plot(self, xmin, xmax, idx_input=0, idx_output=0, points=100,
             **kwargs):
        """Call method |ANN.plot| of all |anntools.ANN| objects
        handled bythe actual |anntools.SeasonalANN| object.
        """
        for toy, ann in self:
            ann.plot(xmin, xmax,
                     idx_input=idx_input, idx_output=idx_output,
                     points=points,
                     label=str(toy))
        pyplot.legend()

    def __getattribute__(self, name):
        if name.startswith('toy_'):
            try:
                try:
                    return self._toy2ann[timetools.TOY(name)]
                except KeyError:
                    raise AttributeError(
                        'No neural network is registered under '
                        'a TOY object named `%s`.'
                        % timetools.TOY(name))
            except BaseException:
                objecttools.augment_excmessage(
                    'While trying to look up for a neural network '
                    'handled by the seasonal neural network collection '
                    '`%s` of element `%s` based on name `%s`'
                    % (self.name, objecttools.devicename(self), name))
        else:
            return object.__getattribute__(self, name)

    def __setattr__(self, name, value):
        if name.startswith('toy_'):
            try:
                if not isinstance(value, abctools.ANNABC):
                    raise TypeError(
                        '%s has been given, but a value of type '
                        '`ANN` is required.'
                        % objecttools.value_of_type(value).capitalize())
                self._toy2ann[timetools.TOY(name)] = value
                self.refresh()
            except BaseException:
                objecttools.augment_excmessage(
                    'While trying to assign a new neural network to '
                    'the seasonal neural network collection `%s` of '
                    'element `%s` based on name `%s`'
                    % (self.name, objecttools.devicename(self), name))
        else:
            object.__setattr__(self, name, value)

    def __delattr__(self, name):
        if name.startswith('toy_'):
            try:
                try:
                    del self._toy2ann[timetools.TOY(name)]
                except KeyError:
                    raise AttributeError(
                        'No neural network is registered under '
                        'a TOY object named `%s`.'
                        % timetools.TOY(name))
                self.refresh()
            except BaseException:
                objecttools.augment_excmessage(
                    'While trying to remove a new neural network from '
                    'the seasonal neural network collection `%s` of '
                    'element `%s` based on name `%s`'
                    % (self.name, objecttools.devicename(self), name))
        else:
            object.__delattr__(self, name)

    def __iter__(self):
        for toy, ann_ in sorted(self._toy2ann.items()):
            yield (toy, ann_)

    def __repr__(self):
        if not self:
            return self.name+'()'
        elif (len(self) == 1) and (self.toys[0] == timetools.TOY('1_1_0_0_0')):
            return self.anns[0].assignrepr('%s(' % self.name) + ')'
        lines = []
        blanks = ' '*(len(self.name)+1)
        for idx, (toy, ann_) in enumerate(self):
            if idx == 0:
                prefix = '%s(%s=' % (self.name, toy)
            else:
                prefix = '%s%s=' % (blanks, toy)
            lines.append(ann_.assignrepr(prefix))
        lines[-1] += ')'
        return ',\n'.join(lines)

    def __len__(self):
        return len(self._toy2ann)

    def __dir__(self):
        """
        >>> from hydpy import SeasonalANN, ann
        >>> seasonalann = SeasonalANN()
        >>> seasonalann(ann(nmb_inputs=1, nmb_neurons=(1,), nmb_outputs=1,
        ...                 weights_input=0.0, weights_output=0.0,
        ...                 intercepts_hidden=0.0, intercepts_output=1.0))
        >>> from hydpy.core.objecttools import assignrepr_values
        >>> print(assignrepr_values(sorted(dir(seasonalann)), '', 70))
        NDIM, SPAN, TIME, TYPE, anns, connect, fastaccess, inputs, name,
        nmb_inputs, nmb_outputs, outputs, parameterstep, plot,
        process_actual_input, ratios, refresh, shape, simulationstep, subpars,
        toy_1_1_0_0_0, toys, verify
        """
        return objecttools.dir_(self) + [str(toy) for toy in self.toys]


abctools.ParameterABC.register(SeasonalANN)
abctools.SeasonalANNABC.register(SeasonalANN)


autodoctools.autodoc_module()

# -*- coding: utf-8 -*-
"""This module provides additional features for module |iuhtools|,
related to Autoregressive-Moving Average (ARMA) models."""
# import...
# ...from standard library
from __future__ import division, print_function
import itertools
import warnings
# ...from site-packages
import numpy
from scipy import integrate
from matplotlib import pyplot
# ...from HydPy
from hydpy import pub
from hydpy.core import autodoctools
from hydpy.core import objecttools
from hydpy.auxs import statstools


class MA(object):
    """Moving Average Model.

    The MA coefficients can be set manually:

    >>> from hydpy import MA
    >>> ma = MA(coefs=(0.8, 0.2))
    >>> ma
    MA(coefs=(0.8, 0.2))
    >>> ma.coefs = 0.2, 0.8
    >>> ma
    MA(coefs=(0.2, 0.8))

    Otherwise they are determined by method |MA.update_coefs|.
    But this requires that a integrable function object is given.
    Usually, this function object is a |IUH| subclass object, but
    (as in the following example) other function objects defining
    instantaneuous unit hydrographs are accepted.  However, they
    should be well-behaved (e.g. be relatively smooth, unimodal,
    strictly positive, unity integral surface in the positive range).

    For educational purposes, some discontinuous functions are applied in
    the following. One can suppress the associated warning messages with
    the following commands:

    >>> import warnings
    >>> from scipy import integrate
    >>> warnings.filterwarnings(
    ... 'ignore', category=integrate.IntegrationWarning)

    The first example is a simple rectangle impuls:

    >>> ma = MA(iuh=lambda x: 0.05 if x < 20.0 else 0.0)
    >>> ma.iuh.moment1 = 10.0
    >>> ma
    MA(coefs=(0.025, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05,
              0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05,
              0.025))

    The number of the coefficients can be modified by changing the
    class attribute |MA.smallest_coeff|:

    >>> ma.smallest_coeff = 0.03
    >>> ma.update_coefs()
    >>> ma
    MA(coefs=(0.025641, 0.051282, 0.051282, 0.051282, 0.051282, 0.051282,
              0.051282, 0.051282, 0.051282, 0.051282, 0.051282, 0.051282,
              0.051282, 0.051282, 0.051282, 0.051282, 0.051282, 0.051282,
              0.051282, 0.051282))


    The first two central moments of the time delay are a usefull measure for
    describing the operation of a MA model:

    >>> ma = MA(iuh=lambda x: 1.0 if x < 1.0 else 0.0)
    >>> ma.iuh.moment1 = 0.5
    >>> ma
    MA(coefs=(0.5, 0.5))
    >>> from hydpy import round_
    >>> round_(ma.moments, 6)
    0.5, 0.5

    The first central moment is the weigthed time delay (mean lag time).
    The second central moment is the weighted mean deviation from the
    mean lag time (diffusion).

    MA objects can return the turning point in the recession part of their
    MA coefficients.  This can be demonstrated for the right side of the
    probability density function of the normal distribution with zero
    mean and a standard deviation (turning point) of 10:

    >>> from scipy import stats
    >>> ma = MA(iuh=lambda x: 2.0*stats.norm.pdf(x, 0.0, 2.0))
    >>> ma.iuh.moment1 = 1.35
    >>> ma
    MA(coefs=(0.195417, 0.346659, 0.24189, 0.13277, 0.057318, 0.019458,
              0.005193, 0.001089, 0.00018, 0.000023, 0.000002, 0.0, 0.0))
    >>> round_(ma.turningpoint)
    2, 0.24189

    Note that the first returned value is the index of the the MA coefficient
    closest to the turning point, and not a high precision estimate of
    the real turning point of the instantaneous unit hydrograph.

    You can also use the following ploting command to verify the position of
    the turning point, which is printed as a red dot.

    >>> ma.plot(threshold=0.9)

    .. testsetup::

        >>> from matplotlib import pyplot
        >>> pyplot.close()

    The turning point detection also works for functions which include
    both a rising and a falling limb.  This can be shown shifting the
    normal distribution to the right:

    >>> ma.iuh = lambda x: 1.02328*stats.norm.pdf(x, 4.0, 2.0)
    >>> ma.iuh.moment1 = 3.94
    >>> ma.update_coefs()
    >>> ma
    MA(coefs=(0.019322, 0.067931, 0.12376, 0.177364, 0.199966, 0.177364,
              0.12376, 0.067931, 0.029326, 0.009956, 0.002657, 0.000557,
              0.000092, 0.000012, 0.000001, 0.0, 0.0))
    >>> round_(ma.turningpoint)
    6, 0.12376

    When no turning point can be detected, an error is raised:

    >>> ma.coefs = 1.0, 1.0, 1.0
    >>> ma.turningpoint
    Traceback (most recent call last):
    ...
    RuntimeError: Not able to detect a turning point in the impulse response \
defined by the MA coefficients 1.0, 1.0, 1.0.

    The next example requires reactivating the warning suppressed above
    (and under Python 2.7 some registry clearing):

    >>> warnings.filterwarnings(
    ...     'error', category=integrate.IntegrationWarning)
    >>> from scipy.integrate import quadpack
    >>> quadpack.__warningregistry__ = {}

    The MA coefficients need to be approximated numerically.  For very
    spiky response function, the underlying integration algorithm might
    fail.  Then it is assumed that the complete mass of the response
    function is placed at a single delay time, defined by the property
    `moment1` of the instantaneous unit hydrograph.  Hopefully, this
    leads to plausible results.  However, we raise an additional warning
    message to allow users to determine the coefficients by a different
    approach :

    >>> ma.iuh = lambda x: 10.0 if 4.2 < x <= 4.3 else 0.0
    >>> ma.iuh.moment1 = 4.25
    >>> ma.update_coefs()
    Traceback (most recent call last):
    ...
    UserWarning: During the determination of the MA coefficients \
corresponding to the instantaneous unit hydrograph ... a numerical \
integration problem occured.  \
Please check the calculated coefficients: 0.0, 0.0, 0.0, 0.0, 0.75, 0.25.
    >>> ma
    MA(coefs=(0.0, 0.0, 0.0, 0.0, 0.75, 0.25))
    """

    smallest_coeff = 1e-9
    """Smalles MA coefficient to be determined at the end of the response."""

    _coefs = None

    def __init__(self, iuh=None, coefs=None):
        self.iuh = iuh
        if coefs is not None:
            self.coefs = coefs

    def _get_coefs(self):
        """|numpy.ndarray| containing all MA coefficents."""
        if self._coefs is None:
            self.update_coefs()
        return self._coefs

    def _set_coefs(self, values):
        self._coefs = numpy.array(values, ndmin=1, dtype=float)

    def _del_coefs(self):
        self._coefs = None

    coefs = property(_get_coefs, _set_coefs, _del_coefs)

    @property
    def order(self):
        """MA order."""
        return len(self.coefs)

    def _quad(self, dt, t):
        return integrate.quad(
            self.iuh, max(t-1.+dt, 0.), t+dt)[0]

    def update_coefs(self):
        """(Re)calculate the MA coefficients based on the instantaneous
        unit hydrograph."""
        coefs = []
        sum_coefs = 0.
        moment1 = self.iuh.moment1
        for t in itertools.count(0., 1.):
            points = (moment1 % 1,) if t <= moment1 <= (t+2.) else ()
            try:
                coef = integrate.quad(
                    self._quad, 0., 1., args=(t,), points=points)[0]
            except integrate.IntegrationWarning:
                idx = int(moment1)
                coefs = numpy.zeros(idx+2, dtype=float)
                weight = (moment1-idx)
                coefs[idx] = (1.-weight)
                coefs[idx+1] = weight
                self.coefs = coefs
                warnings.warn(
                    'During the determination of the MA coefficients '
                    'corresponding to the instantaneous unit hydrograph '
                    '`%s` a numerical integration problem occured.  '
                    'Please check the calculated coefficients: %s.'
                    % (repr(self.iuh), objecttools.repr_values(coefs)))
                break   # pragma: no cover
            sum_coefs += coef
            if (sum_coefs > .9) and (coef < self.smallest_coeff):
                coefs = numpy.array(coefs)
                coefs /= numpy.sum(coefs)
                self.coefs = coefs
                break
            else:
                coefs.append(coef)

    @property
    def turningpoint(self):
        """Turning point (index and value tuple) in the recession part of the
        MA approximation of the instantaneous unit hydrograph."""
        coefs = self.coefs
        old_dc = coefs[1]-coefs[0]
        for idx in range(self.order-2):
            new_dc = coefs[idx+2]-coefs[idx+1]
            if (old_dc < 0.) and (new_dc > old_dc):
                return idx, coefs[idx]
            else:
                old_dc = new_dc
        raise RuntimeError(
            'Not able to detect a turning point in the impulse response '
            'defined by the MA coefficients %s.'
            % objecttools.repr_values(coefs))

    @property
    def delays(self):
        """Time delays related to the individual MA coefficients."""
        return numpy.arange(self.order, dtype=float)

    @property
    def moments(self):
        """The first two time delay weighted statistical moments of the
        MA coefficients."""
        moment1 = statstools.calc_mean_time(self.delays, self.coefs)
        moment2 = statstools.calc_mean_time_deviation(
            self.delays, self.coefs, moment1)
        return numpy.array([moment1, moment2])

    def plot(self, threshold=None, **kwargs):
        """Barplot of the MA coefficients."""
        try:
            # Works under matplotlib 3.
            pyplot.bar(x=self.delays+.5, height=self.coefs,
                       width=1., fill=False, **kwargs)
        except TypeError:   # pragma: no cover
            # Works under matplotlib 2.
            pyplot.bar(left=self.delays+.5, height=self.coefs,
                       width=1., fill=False, **kwargs)
        pyplot.xlabel('time')
        pyplot.ylabel('response')
        if threshold is not None:
            cumsum = numpy.cumsum(self.coefs)
            idx = numpy.where(cumsum > threshold*cumsum[-1])[0][0]
            pyplot.xlim(0., idx)
        idx, value = self.turningpoint
        pyplot.plot(idx, value, 'ro')

    def __repr__(self):
        return objecttools.assignrepr_tuple(self.coefs, 'MA(coefs=', 70) + ')'


class ARMA(object):
    """Autoregressive-Moving Average model.

    All ARMA coefficients can be set manually:

    >>> from hydpy import MA, ARMA
    >>> arma = ARMA(ar_coefs=(0.5,), ma_coefs=(0.3, 0.2))
    >>> arma.coefs
    (array([ 0.5]), array([ 0.3,  0.2]))
    >>> arma
    ARMA(ar_coefs=(0.5,),
         ma_coefs=(0.3, 0.2))

    >>> arma.ar_coefs = ()
    >>> arma.ma_coefs = range(20)
    >>> arma
    ARMA(ar_coefs=(),
         ma_coefs=(0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0,
                   11.0, 12.0, 13.0, 14.0, 15.0, 16.0, 17.0, 18.0, 19.0))

    Otherwise they are determined by method |ARMA.update_coefs|.
    But this requires that a |MA| object is given.  Let us use
    the MA model based on the shifted normal distribution of the
    documentation on class |MA| as an example:

    >>> from scipy import stats
    >>> ma = MA(iuh=lambda x: 1.02328*stats.norm.pdf(x, 4.0, 2.0))
    >>> ma.iuh.moment1 = 3.94
    >>> arma = ARMA(ma_model=ma)
    >>> arma
    ARMA(ar_coefs=(0.680483, -0.228511, 0.047283, -0.006022, 0.000377),
         ma_coefs=(0.019322, 0.054783, 0.08195, 0.107757, 0.104458,
                   0.07637, 0.041095, 0.01581, 0.004132, 0.000663,
                   0.00005))

    To verify that the ARMA model approximates the MA model with sufficient
    accuracy, one can check the central moments of their responses to a
    standard delta time impulse:

    >>> from hydpy import round_
    >>> round_(ma.moments)
    4.110496, 1.926798
    >>> round_(arma.moments)
    4.110496, 1.926798

    On can check the accuray of the approximation directly via the
    property |ARMA.dev_moments|, which is the sum of the absolute values
    of the deviations of both methods:

    >>> round_(arma.dev_moments)
    0.0

    For the first six digits, there is no difference.  However, the total
    number of coefficients is only reduced by one:

    >>> ma.order
    17
    >>> arma.order
    (5, 11)

    To reduce the determined number or AR coefficients, one can set a higher
    AR-related tolerance value:

    >>> arma.max_rel_rmse = 1e-3
    >>> arma.update_coefs()
    >>> arma
    ARMA(ar_coefs=(0.788899, -0.256436, 0.034256),
         ma_coefs=(0.019322, 0.052688, 0.075125, 0.096488, 0.089453,
                   0.060854, 0.029041, 0.008929, 0.001397, 0.000001,
                   -0.000004, 0.00001, -0.000008, -0.000009, -0.000004,
                   -0.000001))

    The number of AR coeffcients is actually reduced.  However, the are
    now even more MA coefficients, possibly trying to compensate the lower
    accuracy of the AR coefficients, and there is a slight decrease in
    the precision of the moments:

    >>> arma.order
    (3, 16)
    >>> round_(arma.moments)
    4.110497, 1.926804
    >>> round_(arma.dev_moments)
    0.000007

    To also reduce the number of MA coefficients, one can set a higher
    MA-related tolerance value:

    >>> arma.max_dev_coefs = 1e-3
    >>> arma.update_coefs()
    >>> arma
    ARMA(ar_coefs=(0.788888, -0.256432, 0.034255),
         ma_coefs=(0.019321, 0.052687, 0.075124, 0.096486, 0.089452,
                   0.060853, 0.02904, 0.008929, 0.001397))

    Now the total number of is in fact decreased, and the loss in accuracy
    is still small:

    >>> arma.order
    (3, 9)
    >>> round_(arma.moments)
    4.110794, 1.927625
    >>> round_(arma.dev_moments)
    0.001125

    Further relaxing the tolerance values results in even less coefficients,
    but also in some slightly negative responses to a standard impulse:

    >>> arma.max_rel_rmse = 1e-2
    >>> arma.max_dev_coefs = 1e-2
    >>> arma.update_coefs()
    Traceback (most recent call last):
    ...
    UserWarning: Note that the smallest response to a standard impulse of the \
determined ARMA model is negative (`-0.000316`).
    >>> arma
    ARMA(ar_coefs=(0.736954, -0.166457),
         ma_coefs=(0.01946, 0.05418, 0.077804, 0.098741, 0.091295,
                   0.060797, 0.027226))
    >>> arma.order
    (2, 7)
    >>> from hydpy import print_values
    >>> print_values(arma.response)
    0.01946, 0.068521, 0.125062, 0.1795, 0.202761, 0.180343, 0.12638,
    0.063117, 0.025477, 0.008269, 0.001853, -0.000011, -0.000316,
    -0.000231, -0.000118, -0.000048, -0.000016

    It seems to be hard to find a parameter efficient approximation to the
    MA model in the given example. Generally, approximating ARMA models to MA
    models is more beneficial when functions with long tails are involved.
    The most extreme example would be a simple exponential decline:

    >>> import numpy
    >>> ma = MA(iuh=lambda x: 0.1*numpy.exp(-0.1*x))
    >>> ma.iuh.moment1 = 6.932
    >>> arma = ARMA(ma_model=ma)

    In the given example a number of 185 MA coefficients can be reduced to a
    total number of three ARMA coefficients with no relevant loss of accuracy:

    >>> ma.order
    185
    >>> arma.order
    (1, 2)
    >>> round_(arma.dev_moments)
    0.0

    Use the following plotting command to see why 2 MA coeffcients instead of
    one are required in the above example:

    >>> arma.plot(threshold=0.9)

    Decreasing the tolerance values too much results in the following errors:

    >>> arma.max_dev_coefs = 0.0
    >>> arma.update_coefs()
    Traceback (most recent call last):
    ...
    RuntimeError: Method `update_ma_coefs` is not able to determine the MA \
coefficients of the ARMA model with the desired accuracy.  You can set the \
tolerance value ´max_dev_coefs` to a higher value.  An accuracy of \
`0.000000000925` has been reached using `185` MA coefficients.
    >>> arma.max_rel_rmse = 0.0
    >>> arma.update_coefs()
    Traceback (most recent call last):
    ...
    RuntimeError: Method `update_ar_coefs` is not able to determine the AR \
coefficients of the ARMA model with the desired accuracy.  You can either \
set the tolerance value `max_rel_rmse` to a higher value or increase the \
allowed `max_ar_order`.  An accuracy of `0.0` has been reached using `10` \
coefficients.
    """

    max_ar_order = 10
    """Maximum number of AR coefficients that are to be determined by method
    |ARMA.update_coefs|."""

    max_rel_rmse = 1e-6
    """Maximum relative root mean squared error to be accepted by method
    |ARMA.update_coefs|."""

    max_dev_coefs = 1e-6
    """Maximum deviation of the sum of all coefficents from one to be accepted
    by method |ARMA.update_coefs|."""

    _ma_coefs = None
    _ar_coefs = None

    def __init__(self, ma_model=None, ar_coefs=None, ma_coefs=None):
        self.ma = ma_model
        if ar_coefs is not None:
            self.ar_coefs = ar_coefs
        if ma_coefs is not None:
            self.ma_coefs = ma_coefs
        self._rel_rmse = None

    @property
    def rel_rmse(self):
        """Relative root mean squared error the last time achieved by method
        |ARMA.update_coefs|."""
        return self._rel_rmse

    def _get_ar_coefs(self):
        """The AR coefficients of the AR model."""
        if self._ar_coefs is None:
            self.update_ar_coefs()
        return self._ar_coefs

    def _set_ar_coefs(self, values):
        self._ar_coefs = numpy.array(values, ndmin=1, dtype=float)

    def _del_ar_coefs(self):
        self._ar_coefs = None

    ar_coefs = property(_get_ar_coefs, _set_ar_coefs, _del_ar_coefs)

    def _get_ma_coefs(self):
        """The MA coefficients of the ARMA model."""
        if self._ma_coefs is None:
            self.update_ma_coefs()
        return self._ma_coefs

    def _set_ma_coefs(self, values):
        self._ma_coefs = numpy.array(values, ndmin=1, dtype=float)

    def _del_ma_coefs(self):
        self._ma_coefs = None

    ma_coefs = property(_get_ma_coefs, _set_ma_coefs, _del_ma_coefs)

    @property
    def coefs(self):
        """Tuple containing both the AR and the MA coefficients."""
        return (self.ar_coefs, self.ma_coefs)

    @property
    def ar_order(self):
        """Number of AR coefficients."""
        return len(self.ar_coefs)

    @property
    def ma_order(self):
        """Number of MA coefficients"""
        return len(self.ma_coefs)

    @property
    def order(self):
        """Number of both the AR and the MA coefficients."""
        return (self.ar_order, self.ma_order)

    def update_coefs(self):
        """Determine both the AR and the MA coefficients."""
        self.update_ar_coefs()
        self.update_ma_coefs()

    @property
    def effective_max_ar_order(self):
        """The maximum number of AR coefficients that shall or can be
        determined.

        It is the minimum of |ARMA.max_ar_order| and the number of
        coefficients of the pure |MA| after their turning point.
        """
        return min(self.max_ar_order, self.ma.order-self.ma.turningpoint[0]-1)

    def update_ar_coefs(self):
        """Determine the AR coefficients.

        The number of AR coefficients is subsequently increased until the
        required precision |ARMA.max_rel_rmse| is reached.  Otherwise,
        a |RuntimeError| is raised.
        """
        del self.ar_coefs
        for ar_order in range(1, self.effective_max_ar_order+1):
            self.calc_all_ar_coefs(ar_order, self.ma)
            if self._rel_rmse < self.max_rel_rmse:
                break
        else:
            with pub.options.reprdigits(12):
                raise RuntimeError(
                    'Method `update_ar_coefs` is not able to determine '
                    'the AR coefficients of the ARMA model with the desired '
                    'accuracy.  You can either set the tolerance value '
                    '`max_rel_rmse` to a higher value or increase the '
                    'allowed `max_ar_order`.  An accuracy of `%s` has been '
                    'reached using `%d` coefficients.'
                    % (objecttools.repr_(self._rel_rmse), ar_order))

    @property
    def dev_moments(self):
        """Sum of the absolute deviations between the central moments of the
        instantaneous unit hydrograph and the ARMA approximation."""
        return numpy.sum(numpy.abs(self.moments-self.ma.moments))

    def norm_coefs(self):
        """Multiply all coefficients by the same factor, so that their sum
        becomes one."""
        sum_coefs = self.sum_coefs
        self.ar_coefs /= sum_coefs
        self.ma_coefs /= sum_coefs

    @property
    def sum_coefs(self):
        """The sum of all AR and MA coefficients"""
        return numpy.sum(self.ar_coefs) + numpy.sum(self.ma_coefs)

    @property
    def dev_coefs(self):
        """Absolute deviation of |ARMA.sum_coefs| from one."""
        return abs(self.sum_coefs-1.)

    def calc_all_ar_coefs(self, ar_order, ma_model):
        """Determine the AR coeffcients based on a least squares approach.

        The argument `ar_order` defines the number of AR coefficients to be
        determined.  The argument `ma_order` defines a pure |MA| model.
        The least squares approach is applied on all those coefficents of the
        pure MA model, which are associated with the part of the recession
        curve behind its turning point.

        The attribute |ARMA.rel_rmse| is updated with the resulting
        relative root mean square error.
        """
        turning_idx, _ = ma_model.turningpoint
        values = ma_model.coefs[turning_idx:]
        self.ar_coefs, residuals = numpy.linalg.lstsq(
            self.get_a(values, ar_order),
            self.get_b(values, ar_order),
            rcond=-1)[:2]
        if len(residuals) == 1:
            self._rel_rmse = numpy.sqrt(residuals[0])/numpy.sum(values)
        else:
            self._rel_rmse = 0.

    @staticmethod
    def get_a(values, n):
        """Extract the independent variables of the given values and return
        them as a matrix with n columns in a form suitable for the least
        squares approach applied in method |ARMA.update_ar_coefs|.
        """
        m = len(values)-n
        a = numpy.empty((m, n), dtype=float)
        for i in range(m):
            i0 = i-1 if i > 0 else None
            i1 = i+n-1
            a[i] = values[i1:i0:-1]
        return numpy.array(a)

    @staticmethod
    def get_b(values, n):
        """Extract the dependent variables of the values in a vector with n
        entries in a form suitable for the least squares approach applied in
        method |ARMA.update_ar_coefs|.
        """
        return numpy.array(values[n:])

    def update_ma_coefs(self):
        """Determine the MA coefficients.

        The number of MA coefficients is subsequently increased until the
        required precision |ARMA.max_dev_coefs| is reached.  Otherwise,
        a |RuntimeError| is raised.
        """
        self.ma_coefs = []
        for ma_order in range(1, self.ma.order+1):
            self.calc_next_ma_coef(ma_order, self.ma)
            if self.dev_coefs < self.max_dev_coefs:
                self.norm_coefs()
                break
        else:
            with pub.options.reprdigits(12):
                raise RuntimeError(
                    'Method `update_ma_coefs` is not able to determine the MA '
                    'coefficients of the ARMA model with the desired accuracy.'
                    '  You can set the tolerance value ´max_dev_coefs` to a '
                    'higher value.  An accuracy of `%s` has been reached '
                    'using `%d` MA coefficients.'
                    % (objecttools.repr_(self.dev_coefs), ma_order))
        if numpy.min(self.response) < 0.:
            warnings.warn(
                'Note that the smallest response to a standard impulse of the '
                'determined ARMA model is negative (`%s`).'
                % objecttools.repr_(numpy.min(self.response)))

    def calc_next_ma_coef(self, ma_order, ma_model):
        """Determine the MA coefficients of the ARMA model based on its
        predetermined AR coefficients and the MA ordinates of the given
        |MA| model.

        The MA coefficients are determined one at a time, beginning with the
        first one.  Each ARMA MA coefficient in set in a manner that allows
        for the exact reproduction of the equivalent pure MA coefficient with
        all relevant ARMA coefficients.
        """
        idx = ma_order-1
        coef = ma_model.coefs[idx]
        for jdx, ar_coef in enumerate(self.ar_coefs):
            zdx = idx-jdx-1
            if zdx >= 0:
                coef -= ar_coef*ma_model.coefs[zdx]
        self.ma_coefs = numpy.concatenate((self.ma_coefs, [coef]))

    @property
    def response(self):
        """Return the response to a standard dt impulse."""
        values = []
        sum_values = 0.
        idx = 0
        ma_coefs = self.ma_coefs
        ar_coefs = self.ar_coefs
        ma_order = self.ma_order
        for idx in range(len(self.ma.delays)):
            value = 0.
            if idx < ma_order:
                value += ma_coefs[idx]
            for jdx, ar_coef in enumerate(ar_coefs):
                zdx = idx-jdx-1
                if zdx >= 0:
                    value += ar_coef*values[zdx]
            values.append(value)
            sum_values += value
        return numpy.array(values)

    @property
    def moments(self):
        """The first two time delay weighted statistical moments of the
        ARMA response."""
        timepoints = self.ma.delays
        response = self.response
        moment1 = statstools.calc_mean_time(timepoints, response)
        moment2 = statstools.calc_mean_time_deviation(
            timepoints, response, moment1)
        return numpy.array([moment1, moment2])

    def plot(self, threshold=None, **kwargs):
        """Barplot of the ARMA response."""
        try:
            # Works under matplotlib 3.
            pyplot.bar(x=self.ma.delays+.5, height=self.response,
                       width=1., fill=False, **kwargs)
        except TypeError:   # pragma: no cover
            # Works under matplotlib 2.
            pyplot.bar(left=self.ma.delays+.5, height=self.response,
                       width=1., fill=False, **kwargs)
        pyplot.xlabel('time')
        pyplot.ylabel('response')
        if threshold is not None:
            cumsum = numpy.cumsum(self.response)
            idx = numpy.where(cumsum > threshold*cumsum[-1])[0][0]
            pyplot.xlim(0., idx)

    def __repr__(self):
        return '%s,\n%s)' % (objecttools.assignrepr_tuple(self.ar_coefs,
                                                          'ARMA(ar_coefs=',
                                                          70),
                             objecttools.assignrepr_tuple(self.ma_coefs,
                                                          '     ma_coefs=',
                                                          70))


autodoctools.autodoc_module()

# -*- coding: utf-8 -*-
"""This module supports modelling based on instantaneous unit hydrographs.

This module implements some abstract descriptor classes, metaclasses and base
classes.  If you are just interested in applying a certain instantaneous
unit hydrograph (iuh) function or if you want to implement an additional
iuh, see the examples or the source code of class
|TranslationDiffusionEquation|.
"""
# import...
# ...from standard library
from __future__ import division, print_function
import itertools
# ...from site-packages
import numpy
from scipy import special
from matplotlib import pyplot
# ...from Hydpy
from hydpy.core import autodoctools
from hydpy.core import objecttools
from hydpy.auxs import statstools
from hydpy.auxs import armatools


class IUH_Parameter(object):
    """Descriptor base class for |PrimaryParameter| and |SecondaryParameter|.

    The first initialization argument is the parameters name.  Optionally,
    an alternative type (the default type is |float|) and a documentation
    string can be passed.
    """

    def __init__(self, name, type_=float, doc=None):
        self.name = name
        self._name = '_'+name
        self.type_ = type_
        self.__doc__ = ('Instantaneous unit hydrograph parameter %s.' % name
                        if doc is None else str(doc))

    def __get__(self, obj, type_=None):
        return self if obj is None else getattr(obj, self._name, None)

    def _convert_type(self, value):
        try:
            return self.type_(value)
        except BaseException:
            raise TypeError(
                'The value `%s` of type `%s` could not be converted to type '
                '`%s` of the instantaneous unit hydrograph parameter `%s`.'
                % (value, objecttools.classname(value),
                   objecttools.classname(self.type_), self.name))


class PrimaryParameter(IUH_Parameter):
    """Descriptor base class for parameters of instantaneous unit hydrograph
    functions to be defined by the user.

    When a primary parameter value is set or deleted, the master instance is
    instructed to |IUH.update| all secondary parameter values.
    """

    def __set__(self, obj, value):
        value = self._convert_type(value)
        setattr(obj, self._name, value)
        obj.update()

    def __delete__(self, obj):
        setattr(obj, self._name, None)
        obj.update()


class SecondaryParameter(IUH_Parameter):
    """Descriptor base class for parameters of instantaneous unit hydrograph
    functions which can be determined automatically."""

    def __set__(self, obj, value):
        value = self._convert_type(value)
        setattr(obj, self._name, value)

    def __delete__(self, obj):
        setattr(obj, self._name, None)


class MetaIUH(type):
    """Metaclass for class |IUH|.

    For storing |PrimaryParameter| and |SecondaryParameter| in separate
    dictionaries.
    """

    def __new__(cls, name, parents, dict_):
        primary_parameters = {}
        secondary_parameters = {}
        for key, value in dict_.items():
            if isinstance(value, PrimaryParameter):
                primary_parameters[key] = value
            elif isinstance(value, SecondaryParameter):
                secondary_parameters[key] = value
        dict_['_primary_parameters'] = primary_parameters
        dict_['_secondary_parameters'] = secondary_parameters
        return type.__new__(cls, name, parents, dict_)


# Just for making MetaIUH the type of class IUH both in Python 2 and 3:
_MetaIUH = MetaIUH('_MetaIUH', (), {})


class IUH(_MetaIUH):
    """Base class for instantaneous unit hydrograph function objects.

    See class |TranslationDiffusionEquation| for explanations and
    application examples.

    For developers: The string representation does also work for
    parameter-free |IUH| subclasses:

    >>> from hydpy.auxs.iuhtools import IUH
    >>> class Test(IUH):
    ...     pass
    >>> Test()
    Test()
    """

    dt_response = 1e-2
    """Relative stepsize for plotting and analyzing iuh functions."""

    smallest_response = 1e-9
    """Smallest value taken into account for plotting and analyzing iuh
    functions."""

    def __init__(self, **kwargs):
        self.ma = armatools.MA(self)
        self.arma = armatools.ARMA(ma_model=self.ma)
        if kwargs:
            self.set_primary_parameters(**kwargs)

    def set_primary_parameters(self, **kwargs):
        """Set all primary parameters at once."""
        given = sorted(kwargs.keys())
        required = sorted(self._primary_parameters)
        if given == required:
            for (key, value) in kwargs.items():
                setattr(self, key, value)
        else:
            raise ValueError(
                'When passing primary parameter values as initialization '
                'arguments of the instantaneous unit hydrograph class `%s`, '
                'or when using method `set_primary_parameters, one has to '
                'to define all values at once via keyword arguments.  '
                'But instead of the primary parameter names `%s` the '
                'following keywords were given: %s.'
                % (objecttools.classname(self),
                   ', '.join(required), ', '.join(given)))

    @property
    def are_primary_parameters_complete(self):
        """True/False flag that indicates wheter the values of all primary
        parameters are defined or not."""
        for primpar in self._primary_parameters.values():
            if primpar.__get__(self) is None:
                return False
        return True

    def update(self):
        """Delete the coefficients of the pure MA model and also all MA and
        AR coefficients of the ARMA model.  Also calculate or delete the values
        of all secondary iuh parameters, depending on the completeness of the
        values of the primary parameters.
        """
        del self.ma.coefs
        del self.arma.ma_coefs
        del self.arma.ar_coefs
        if self.are_primary_parameters_complete:
            self.calc_secondary_parameters()
        else:
            for secpar in self._secondary_parameters.values():
                secpar.__delete__(self)

    @property
    def delay_response_series(self):
        """A tuple of two numpy arrays, which hold the time delays and the
        associated iuh values respectively."""
        delays = []
        responses = []
        sum_responses = 0.
        for t in itertools.count(self.dt_response/2., self.dt_response):
            delays.append(t)
            response = self(t)
            responses.append(response)
            sum_responses += self.dt_response*response
            if (sum_responses > .9) and (response < self.smallest_response):
                break
        return numpy.array(delays), numpy.array(responses)

    def plot(self, threshold=None, **kwargs):
        """Plot the instanteneous unit hydrograph.

        The optional argument allows for defining a threshold of the cumulative
        sum uf the hydrograph, used to adjust the largest value of the x-axis.
        It must be a value between zero and one.
        """
        delays, responses = self.delay_response_series
        pyplot.plot(delays, responses, **kwargs)
        pyplot.xlabel('time')
        pyplot.ylabel('response')
        if threshold is not None:
            threshold = numpy.clip(threshold, 0., 1.)
            cumsum = numpy.cumsum(responses)
            idx = numpy.where(cumsum >= threshold*cumsum[-1])[0][0]
            pyplot.xlim(0., delays[idx])

    @property
    def moment1(self):
        """The first time delay weighted statistical moment of the
        instantaneous unit hydrograph."""
        delays, response = self.delay_response_series
        return statstools.calc_mean_time(delays, response)

    @property
    def moment2(self):
        """The second time delay weighted statistical momens of the
        instantaneous unit hydrograph."""
        moment1 = self.moment1
        delays, response = self.delay_response_series
        return statstools.calc_mean_time_deviation(
            delays, response, moment1)

    @property
    def moments(self):
        """The first two time delay weighted statistical moments of the
        instantaneous unit hydrograph."""
        return numpy.array([self.moment1, self.moment2])

    def __repr__(self):
        parts = [objecttools.classname(self), '(']
        for (name, primpar) in sorted(self._primary_parameters.items()):
            value = primpar.__get__(self)
            if value is not None:
                parts.extend([name, '=', objecttools.repr_(value), ', '])
        if parts[-1] == ', ':
            parts[-1] = ')'
        else:
            parts.append(')')
        return ''.join(parts)


class TranslationDiffusionEquation(IUH):
    """An instantaneous unit hydrograph based on the `translation diffusion
    equation`.

    The equation used is a linear approximation of the Saint-Venant
    equations for channel routing:

      :math:`h(t) = \\frac{a}{t \\cdot \\sqrt{\\pi \\cdot t}} \\cdot
      e^{-t \\cdot (a/t-b)^2}`

    with:
      :math:`a = \\frac{x}{2 \\cdot \\sqrt{d}}`

      :math:`b = \\frac{u}{2 \\cdot \\sqrt{d}}`

    There are three primary parameter whichs values need to be defined by
    the user:

    >>> from hydpy import TranslationDiffusionEquation
    >>> tde = TranslationDiffusionEquation(u=5., d=15., x=50.)
    >>> tde
    TranslationDiffusionEquation(d=15.0, u=5.0, x=50.0)

    The values of both secondary parameters are determined automatically:

    >>> from hydpy import round_
    >>> round_((tde.a, tde.b))
    6.454972, 0.645497

    The function can principally be evaluated for time delays larger zero,
    but not for zero time delay, which can cause trouble when applying
    numerical integration algorithms.  This is why we clip the given time
    delay to minimum value of 1e-10 internally.  In most cases (like the
    following), the returned result should be workable for integration
    algorithms:

    >>> import numpy
    >>> round_(tde([0.0, 5.0, 10.0, 15.0, 20.0]))
    0.0, 0.040559, 0.115165, 0.031303, 0.00507

    The first delay weighted central moment of the translation diffusion
    equation corresponds to the time lag (`x`/`u`), the second one to
    wave diffusion:

    >>> round_(tde.moments)
    10.0, 3.464101

    Class |TranslationDiffusionEquation| implements its own property
    `moment1` (used in the example above), which is computationally
    more efficient and robust than the one of its base class |IUH|.
    But both normally, both should return very similar values:

    >>> from hydpy.auxs.iuhtools import IUH
    >>> round_(IUH.moment1.fget(tde))
    10.0


    You can also plot the graph corresponding to the actual parameterization:

    >>> tde.plot(threshold=0.9)

    .. testsetup::

        >>> from matplotlib import pyplot
        >>> pyplot.close()

    All instances of the subclasses of |IUH| provide a pure Moving Average
    and an Autoregressive-Moving Average approximation to the dt standard
    impulse of the instantaneous unit hydrograph function.  In the given
    example, the MA approximation involves 57 coefficients, and the ARMA
    approximation invoves 17 coefficients:

    >>> tde.ma.order
    57
    >>> tde.arma.order
    (3, 14)

    The diffusion of the MA model deviates from the iuh function due to
    aggregation.  For the ARMA model, there is also a slight deviation in
    time delay, as the ARMA model itself is only a approximation of the
    MA model:

    >>> round_(tde.ma.moments)
    10.0, 3.488074
    >>> round_(tde.arma.moments)
    10.000091, 3.488377

    For further information on using MA and ARMA models, read the
    documentation on module |armatools|.

    Changing a primary parameter results in an updating of the secondary
    parameters as well as the MA and the ARMA model:

    >>> tde.x = 5.
    >>> round_((tde.a, tde.b))
    0.645497, 0.645497
    >>> tde.ma.order
    37
    >>> tde.arma.order
    (4, 5)

    As long as the primary parameter values are incomplete, no secondary
    parameter values are available:

    >>> del tde.x
    >>> round_((tde.a, tde.b))
    None, None

    Suitable type conversions are performed when new parameter values are set:

    >>> tde.x = '1.'
    >>> tde.x
    1.0

    It a new value cannot be converted, an error is raised:

    >>> tde.x = 'a'
    Traceback (most recent call last):
    ...
    TypeError: The value `a` of type `str` could not be converted to type \
`float` of the instantaneous unit hydrograph parameter `x`.

    When passing parameter values as initialization arguments or when using
    method `set_primary_parameters`, tests for completeness are performed:

    >>> TranslationDiffusionEquation(u=5.0, d=15.0)
    Traceback (most recent call last):
    ...
    ValueError: When passing primary parameter values as initialization \
arguments of the instantaneous unit hydrograph class \
`TranslationDiffusionEquation`, or when using method \
`set_primary_parameters, one has to to define all values at once via \
keyword arguments.  But instead of the primary parameter names `d, u, x` \
the following keywords were given: d, u.
    """
    u = PrimaryParameter('u', doc='Wave velocity.')
    d = PrimaryParameter('d', doc='Diffusion coefficient.')
    x = PrimaryParameter('x', doc='Routing distance.')
    a = SecondaryParameter('a', doc='Distance related coefficient.')
    b = SecondaryParameter('b', doc='Velocity related coefficient.')

    def calc_secondary_parameters(self):
        """Determine the values of the secondary parameters `a` and `b`."""
        self.a = self.x/(2.*self.d**.5)
        self.b = self.u/(2.*self.d**.5)

    def __call__(self, t):
        t = numpy.array(t)
        t = numpy.clip(t, 1e-10, numpy.inf)
        return self.a/(t*(numpy.pi*t)**.5)*numpy.exp(-t*(self.a/t-self.b)**2)

    @property
    def moment1(self):
        """The first time delay weighted statistical moment of the
        translation diffusion equation."""
        return self.x/self.u


class LinearStorageCascade(IUH):
    """An instantaneous unit hydrograph based on the `linear storage cascade`.

    The equation involves the gamma function, allowing for a fractional number
    of storages:

      :math:`h(t) = c \\cdot (t/k)^{n-1} \\cdot e^{-t/k}`

    with:
      :math:`c = \\frac{1}{k \\cdot \\gamma(n)}`

    After defining the values of the two primary parameters, the function
    object can be applied:

    >>> from hydpy import LinearStorageCascade
    >>> lsc = LinearStorageCascade(n=2.5, k=2.)
    >>> from hydpy import round_
    >>> round_(lsc.c)
    0.376126
    >>> import numpy
    >>> round_(lsc(numpy.array([5., 10., 15., 20.])))
    0.122042, 0.028335, 0.004273, 0.00054

    """
    n = PrimaryParameter('n', doc='Number of linear storages.')
    k = PrimaryParameter(
        'k', doc='Time of concentration of each individual storage.')
    c = SecondaryParameter('c', doc='Proportionality factor.')

    def calc_secondary_parameters(self):
        """Determine the value of the secondary parameter `c`."""
        self.c = 1./(self.k*special.gamma(self.n))

    def __call__(self, t):
        return self.c*(t/self.k)**(self.n-1)*numpy.exp(-t/self.k)

    @property
    def moment1(self):
        """The first time delay weighted statistical moment of the
        linear storage cascade."""
        return self.k*self.n


autodoctools.autodoc_module()

# -*- coding: utf-8 -*-
"""This module provides features for preparing HydPy networks based on
different data."""
# import...
# ...from standard library
from __future__ import division, print_function
# ...from HydPy
from hydpy.core import objecttools
from hydpy.core import devicetools
from hydpy.core import selectiontools


class RiverBasinNumber(str):
    """A single river basin number (Gewässerkennzahl) based on a guideline
    of the German organisation LAWA.

    See `Länderarbeitsgemeinschaft Wasser: Richtlinie für die
    Gebietsbezeichung und die Verschlüsselung von Fließgewässern (1970)`

    Note that zero fill numbers are ignored:

    >>> RiverBasinNumber(1230)
    RiverBasinNumber(123)
    >>> RiverBasinNumber('0123')
    RiverBasinNumber(123)

    Numbers that cannot be interpreted as a river basin numbers, result
    in the following error message:

    >>> RiverBasinNumber('123A')
    Traceback (most recent call last):
    ...
    ValueError: The given value `123A` could not be interpreted as a river \
basin number.
    """

    def __new__(cls, value):
        try:
            return str.__new__(cls, str(int(value)).strip('0'))
        except ValueError:
            raise ValueError('The given value `%s` could not be interpreted '
                             'as a river basin number.' % value)

    @property
    def is_rivermouth(self):
        """Only true if the river basin number ends with `9`.

        >>> from hydpy import RiverBasinNumber
        >>> RiverBasinNumber(129).is_rivermouth
        True
        >>> RiverBasinNumber(123).is_rivermouth
        False
        """
        return self[-1] == '9'

    @property
    def is_mainchannel(self):
        """Only true if the river basin number ends with an odd number.

        >>> from hydpy import RiverBasinNumber
        >>> for number in range(120, 130):
        ...     rbn = RiverBasinNumber(number)
        ...     print(number, rbn, rbn.is_mainchannel)
        120 12 False
        121 121 True
        122 122 False
        123 123 True
        124 124 False
        125 125 True
        126 126 False
        127 127 True
        128 128 False
        129 129 True
        """
        return int(self[-1]) % 2 == 1

    @property
    def is_tributary(self):
        """Only true if the river basin number ends with an even number.

        >>> from hydpy import RiverBasinNumber
        >>> for number in range(120, 130):
        ...     rbn = RiverBasinNumber(number)
        ...     print(number, rbn, rbn.is_tributary)
        120 12 True
        121 121 False
        122 122 True
        123 123 False
        124 124 True
        125 125 False
        126 126 True
        127 127 False
        128 128 True
        129 129 False
        """
        return int(self[-1]) % 2 == 0

    @property
    def _possible_next_last_initial_digits(self):
        return (digit for digit in ('3', '5', '7', '9') if
                digit > self.rstrip('9')[-1])

    @property
    def possible_next_initial_digits(self):
        """Return a tuple of all potential candidates for the next downstream
        river basin number.

        I think, only the first and the last returned candidate should be of
        relevance.  But to return all possible intermediate candidates might
        be safer river basin number specifications I am not taking into
        account at the moment.

        The candidate numbers might be incomplete.  For example, the next
        number downstream of `123` could be `1251` instead of `125`.

        >>> from hydpy import RiverBasinNumber
        >>> for number in range(120, 130):
        ...     rbn = RiverBasinNumber(number)
        ...     print(number, rbn, rbn.possible_next_initial_digits)
        120 12 ('13', '15', '17', '19')
        121 121 ('123', '125', '127', '129')
        122 122 ('123', '125', '127', '129')
        123 123 ('125', '127', '129')
        124 124 ('125', '127', '129')
        125 125 ('127', '129')
        126 126 ('127', '129')
        127 127 ('129',)
        128 128 ('129',)
        129 129 ('13', '15', '17', '19')
        """
        return tuple(self.rstrip('9')[:-1]+digit for digit
                     in self._possible_next_last_initial_digits)

    @property
    def nmb_digits(self):
        """Number of (significant) digits of a river basin number.

        >>> from hydpy import RiverBasinNumber
        >>> RiverBasinNumber(120).nmb_digits
        2
        """
        return len(self)

    def __repr__(self):
        return 'RiverBasinNumber(%s)' % str.__repr__(self)[1:-1]


class RiverBasinNumbers(tuple):
    """A sorted collection of |RiverBasinNumber| objects.

        >>> from hydpy import RiverBasinNumbers
        >>> RiverBasinNumbers((111, 113, 1129, 11269, 1125, 11261,
        ...                    11262, 1123, 1124, 1122, 1121))
        RiverBasinNumbers((111, 1121, 1122, 1123, 1124, 1125, 11261,
                           11262, 11269, 1129, 113))
    """
    def __new__(cls, values):
        _values = tuple(RiverBasinNumber(value) for value in values)
        return tuple.__new__(cls, sorted(_values))

    def _get_next_numbers(self, riverbasinnumber):
        riverbasinnumber = RiverBasinNumber(riverbasinnumber)
        for pdn1 in riverbasinnumber.possible_next_initial_digits:
            neighbours = [rbn for rbn in self if rbn.startswith(pdn1)]
            if neighbours:
                return min(neighbours)
        else:
            return None

    @property
    def next_numbers(self):
        """A tuple of the next downstream river basin numbers.

        The order of the returned numbers corresponds to the order of the
        numbers contained by the |RiverBasinNumbers| object.

        The number of the subcatchment immediately downstream of the outlet
        subcatchment is not known.  The tuple contains a |None| object
        instead (or multiple |None| objects in case of multiple outlets).

        Eventually, not all possible combinations of river basin numbers
        are covered.  Please keep us informed if you notices a problem
        when applying this algorithm on your data.  At least, the algorithm
        works properly on the following test case provided by Michael
        Wagner (TU Dresden):

        .. image:: LAWA_river-basin-bumbers.png

        At first, only the black arrows are considered, exemplifying the basic
        definition of river basin numbers:

        >>> from hydpy import RiverBasinNumbers
        >>> rbns = RiverBasinNumbers((111, 113, 1129, 11269, 1125, 11261,
        ...                           11262, 1123, 1124, 1122, 1121))
        >>> for this_rbn, next_rbn in zip(rbns, rbns.next_numbers):
        ...    print(this_rbn.ljust(6), next_rbn)
        111    113
        1121   1123
        1122   1123
        1123   1125
        1124   1125
        1125   1129
        11261  11269
        11262  11269
        11269  1129
        1129   113
        113    None

        The coloured arrows exemplify the situation, where some additional
        subdivisions become necessary:

        >>> from hydpy import RiverBasinNumbers
        >>> rbns = RiverBasinNumbers((1111, 1119, 113, 1129, 1127,
        ...                           11269, 1125, 11261, 11262, 11239,
        ...                           11231, 1124, 1122, 1121))
        >>> for this_rbn, next_rbn in zip(rbns, rbns.next_numbers):
        ...    print(this_rbn.ljust(6), next_rbn)
        1111   1119
        1119   113
        1121   11231
        1122   11231
        11231  11239
        11239  1125
        1124   1125
        1125   1127
        11261  11269
        11262  11269
        11269  1127
        1127   1129
        1129   113
        113    None
        """
        return tuple(self._get_next_numbers(rbn) for rbn in self)

    def __repr__(self):
        return objecttools.assignrepr_tuple(self,
                                            'RiverBasinNumbers(',
                                            60) + ')'


class RiverBasinNumbers2Selection(object):
    """Class for defining a |Selection| object (consisting of connected
    nodes and elements) based on given |RiverBasinNumber| objects.

    Note that this class is not intended to cover all possible HydPy
    networks.  So it might be necessary to make some adjustments on
    the returned selection, e.g. to define special names for specific
    elements or nodes.

    All examples of the methods and propertys of
    |RiverBasinNumbers2Selection| are based on the river basin numbers
    defined in the documentation on class |RiverBasinNumbers|.
    """

    def __init__(self, numbers):
        self.supplier_prefix = 'land_'
        self.router_prefix = 'stream_'
        self.node_prefix = 'node_'
        self.last_node = 'node_outlet'
        self.selection_name = 'complete'
        rbns = RiverBasinNumbers(numbers)
        self._up2down = dict(tuple_ for tuple_ in zip(rbns, rbns.next_numbers))

    @property
    def _supplier_numbers(self):
        """A tuple of the numbers of all "supplying" basins."""
        return tuple(self._up2down.keys())

    @property
    def _router_numbers(self):
        """A tuple of the numbers of all "routing" basins."""
        return tuple(up for up in self._up2down.keys()
                     if up in self._up2down.values())

    def _get_nodename(self, string):
        return self.node_prefix + string

    def _get_suppliername(self, string):
        return self.supplier_prefix + string

    def _get_routername(self, string):
        return self.router_prefix + string

    @property
    def supplier_elements(self):
        """A |Elements| collection of all "supplying" basins.

        (All river basins are assumed to supply something to the downstream
        basin.)

        >>> from hydpy import RiverBasinNumbers2Selection
        >>> rbns2s = RiverBasinNumbers2Selection(
        ...                            (111, 113, 1129, 11269, 1125, 11261,
        ...                             11262, 1123, 1124, 1122, 1121))

        The following elements are properly connected to the required
        outlet nodes already:

        >>> for element in rbns2s.supplier_elements:
        ...     print(repr(element))
        Element("land_111",
                outlets="node_113")
        Element("land_1121",
                outlets="node_1123")
        Element("land_1122",
                outlets="node_1123")
        Element("land_1123",
                outlets="node_1125")
        Element("land_1124",
                outlets="node_1125")
        Element("land_1125",
                outlets="node_1129")
        Element("land_11261",
                outlets="node_11269")
        Element("land_11262",
                outlets="node_11269")
        Element("land_11269",
                outlets="node_1129")
        Element("land_1129",
                outlets="node_113")
        Element("land_113",
                outlets="node_outlet")

        It is both possible to change the prefix names of the elements
        and nodes, as long as it results in a valid variable name (e.g.
        does not start with a number):

        >>> rbns2s.supplier_prefix = 'a_'
        >>> rbns2s.node_prefix = 'b_'
        >>> rbns2s.supplier_elements
        Elements("a_111", "a_1121", "a_1122", "a_1123", "a_1124", "a_1125",
                 "a_11261", "a_11262", "a_11269", "a_1129", "a_113")
        """
        elements = devicetools.Elements()
        for supplier in self._supplier_numbers:
            element = self._get_suppliername(supplier)
            try:
                outlet = self._get_nodename(self._up2down[supplier])
            except TypeError:
                outlet = self.last_node
            elements += devicetools.Element(element, outlets=outlet)
        return elements

    @property
    def router_elements(self):
        """A |Elements| collection of all "routing" basins.

        (Only river basins with a upstream basin are assumed to route
        something to the downstream basin.)

        >>> from hydpy import RiverBasinNumbers2Selection
        >>> rbns2s = RiverBasinNumbers2Selection(
        ...                            (111, 113, 1129, 11269, 1125, 11261,
        ...                             11262, 1123, 1124, 1122, 1121))

        The following elements are properly connected to the required
        inlet and outlet nodes already:

        >>> for element in rbns2s.router_elements:
        ...     print(repr(element))
        Element("stream_1123",
                inlets="node_1123",
                outlets="node_1125")
        Element("stream_1125",
                inlets="node_1125",
                outlets="node_1129")
        Element("stream_11269",
                inlets="node_11269",
                outlets="node_1129")
        Element("stream_1129",
                inlets="node_1129",
                outlets="node_113")
        Element("stream_113",
                inlets="node_113",
                outlets="node_outlet")

        It is both possible to change the prefix names of the elements
        and nodes, as long as it results in a valid variable name (e.g.
        does not start with a number):

        >>> rbns2s.router_prefix = 'c_'
        >>> rbns2s.node_prefix = 'd_'
        >>> rbns2s.router_elements
        Elements("c_1123", "c_1125", "c_11269", "c_1129", "c_113")
        """
        elements = devicetools.Elements()
        for router in self._router_numbers:
            element = self._get_routername(router)
            inlet = self._get_nodename(router)
            try:
                outlet = self._get_nodename(self._up2down[router])
            except TypeError:
                outlet = self.last_node
            elements += devicetools.Element(
                            element, inlets=inlet, outlets=outlet)
        return elements

    @property
    def elements(self):
        """Both the "supplying" and the "routing" elements."""
        return self.supplier_elements + self.router_elements

    @property
    def nodes(self):
        """A |Nodes| collection of all required nodes.

        >>> from hydpy import RiverBasinNumbers2Selection
        >>> rbns2s = RiverBasinNumbers2Selection(
        ...                            (111, 113, 1129, 11269, 1125, 11261,
        ...                             11262, 1123, 1124, 1122, 1121))

        Note that the required outlet node is added:

        >>> rbns2s.nodes
        Nodes("node_1123", "node_1125", "node_11269", "node_1129", "node_113",
              "node_outlet")

        It is both possible to change the prefix names of the nodes and
        the name of the outlet node separately:

        >>> rbns2s.node_prefix = 'b_'
        >>> rbns2s.last_node = 'l_node'
        >>> rbns2s.nodes
        Nodes("b_1123", "b_1125", "b_11269", "b_1129", "b_113", "l_node")
        """
        return (devicetools.Nodes(self.node_prefix+routers for
                                  routers in self._router_numbers) +
                devicetools.Node(self.last_node))

    @property
    def selection(self):
        """A complete |Selection| object of all "supplying" and "routing"
        elements and required nodes.

        >>> from hydpy import RiverBasinNumbers2Selection
        >>> rbns2s = RiverBasinNumbers2Selection(
        ...                            (111, 113, 1129, 11269, 1125, 11261,
        ...                             11262, 1123, 1124, 1122, 1121))
        >>> rbns2s.selection
        Selection("complete",
                  elements=("land_111", "land_1121", "land_1122", "land_1123",
                            "land_1124", "land_1125", "land_11261",
                            "land_11262", "land_11269", "land_1129",
                            "land_113", "stream_1123", "stream_1125",
                            "stream_11269", "stream_1129", "stream_113"),
                  nodes=("node_1123", "node_1125", "node_11269", "node_1129",
                         "node_113", "node_outlet"))

        Besides the possible modifications on the names of the different
        nodes and elements, the name of the selection can be set differently:

        >>> rbns2s.selection_name = 'sel'
        >>> from hydpy import pub
        >>> with pub.options.ellipsis(1):
        ...     print(repr(rbns2s.selection))
        Selection("sel",
                  elements=("land_111", ...,"stream_113"),
                  nodes=("node_1123", ...,"node_outlet"))
        """
        return selectiontools.Selection(
                            self.selection_name, self.nodes, self.elements)

# -*- coding: utf-8 -*-
"""This module implements features which help to regularize discontinuous
process equations.

.. _Tyralla (2016): http://www.hydrology.ruhr-uni-bochum.de/hydrolgy/mam/download/schriftenreihe_29.pdf

Many hydrological models rely heavily on discontinous equations describing
hydrological processes.  The related "if-else" blocks are often not
theoretically motivated.  Instead, they are thought to ease implementing
ad hoc solutions of different (parts of) process equations without taking
care of the total set of process equations.

There are some reasons to ground new model concepts on mainly continuous
process descriptions. See e.g. `Tyralla (2016)`_ for more a more exhaustive
discussion of this topic.  Nevertheless, one might often want -- least as
a starting point -- to pick single discontinuous but well-established
equations of old model concepts for a new model concept.  The tools
provided by this module can be used to regularize the discontinuities of
such equations.  More concrete, the tools are thought for replacing
discontinous process equations by continuous approximations.

Some of the implemented features are to be applied during model simulations
and are in some other way performance-critical.  These features are defined
in the Cython extension module |smoothutils|.
"""

# import...
# ...from standard library
from __future__ import division, print_function
# ...from site-packages
import numpy
from scipy import optimize
# ...from HydPy
from hydpy.cythons import smoothutils
from hydpy.core import autodoctools


def calc_smoothpar_logistic1(metapar):
    """Return the smoothing parameter corresponding to the given meta
    parameter when using |smooth_logistic1|.

    Calculate the smoothing parameter value corresponding the meta parameter
    value 2.5:

    >>> from hydpy.auxs.smoothtools import calc_smoothpar_logistic1
    >>> smoothpar = calc_smoothpar_logistic1(2.5)

    Using this smoothing parameter value, the output of function
    |smooth_logistic1| differs by 1 % from the related `true`
    discontinuous step function for the input values -2.5 and 2.5
    (which are located at a distance of 2.5 from the position of
    the discontinuity):

    >>> from hydpy.cythons import smoothutils
    >>> from hydpy import round_
    >>> round_(smoothutils.smooth_logistic1(-2.5, smoothpar))
    0.01
    >>> round_(smoothutils.smooth_logistic1(2.5, smoothpar))
    0.99

    For zero or negative meta parameter values, a zero smoothing parameter
    value is returned:

    >>> round_(calc_smoothpar_logistic1(0.0))
    0.0
    >>> round_(calc_smoothpar_logistic1(-1.0))
    0.0
    """
    return max(metapar/numpy.log(99.), 0.)


def _error_smoothpar_logistic2(par, metapar):
    return smoothutils.smooth_logistic2(-metapar, par) - .01


def _smooth_logistic2_derivative(par, metapar):
    return smoothutils.smooth_logistic2_derivative(metapar, par)


def calc_smoothpar_logistic2(metapar):
    """Return the smoothing parameter corresponding to the given meta
    parameter when using |smooth_logistic2|.

    Calculate the smoothing parameter value corresponding the meta parameter
    value 2.5:

    >>> from hydpy.auxs.smoothtools import calc_smoothpar_logistic2
    >>> smoothpar = calc_smoothpar_logistic2(2.5)

    Using this smoothing parameter value, the output of function
    |smooth_logistic2| differs by
    1 % from the related `true` discontinuous step function for the
    input values -2.5 and 2.5 (which are located at a distance of 2.5
    from the position of the discontinuity):

    >>> from hydpy.cythons import smoothutils
    >>> from hydpy import round_
    >>> round_(smoothutils.smooth_logistic2(-2.5, smoothpar))
    0.01
    >>> round_(smoothutils.smooth_logistic2(2.5, smoothpar))
    2.51

    For zero or negative meta parameter values, a zero smoothing parameter
    value is returned:

    >>> round_(calc_smoothpar_logistic2(0.0))
    0.0
    >>> round_(calc_smoothpar_logistic2(-1.0))
    0.0
    """
    if metapar <= 0.:
        return 0.
    else:
        return optimize.newton(_error_smoothpar_logistic2,
                               .3 * metapar**.84,
                               _smooth_logistic2_derivative,
                               args=(metapar,))


def calc_smoothpar_logistic3(metapar):
    """Return the smoothing parameter corresponding to the given meta
    parameter when using |smooth_logistic3|.

    |smooth_logistic3| is only an alias for |smooth_logistic2|.

    Calculate the smoothing parameter value corresponding the meta parameter
    value 2.5:

    >>> from hydpy.auxs.smoothtools import calc_smoothpar_logistic3
    >>> smoothpar = calc_smoothpar_logistic3(2.5)

    Using this smoothing parameter value, the output of function
    |smooth_logistic3| would ideally differs by 1 % from the related
    `true` discontinuous step function for the input values -2.5 and 3.5
    (which are located at a distance of 2.5 from the position of the
    nearest discontinuity):

    >>> from hydpy.cythons import smoothutils
    >>> from hydpy import round_
    >>> round_(smoothutils.smooth_logistic3(-2.5, smoothpar))
    0.009876
    >>> round_(smoothutils.smooth_logistic3(3.5, smoothpar))
    0.990124

    In contrast to the examples shown for functions |smooth_logistic1| and
    |smooth_logistic2|, the smoothing parameter determined for function
    |smooth_logistic3| is not in perfect agreement with the given meta
    parameter.  For most purposes, the resulting error is negligible.
    If one needs a higher accuracy, some iterative refinement should be
    implemented.
    """
    return calc_smoothpar_logistic2(metapar)


def calc_smoothpar_max1(metapar):
    """Return the smoothing parameter corresponding to the given meta
    parameter when using |smooth_max1|.

    |smooth_max1| is only an alias for |smooth_logistic2|.

    Calculate the smoothing parameter value corresponding the meta parameter
    value 2.5:

    >>> from hydpy.auxs.smoothtools import calc_smoothpar_max1
    >>> smoothpar = calc_smoothpar_max1(2.5)

    Using this smoothing parameter value, the output of function
    |smooth_max1| is 0.01 above the usual discontinuous maximum
    function result, if the absolute value of the difference
    between the x and the y value is 2.5:

    >>> from hydpy.cythons import smoothutils
    >>> from hydpy import round_
    >>> round_(smoothutils.smooth_max1(4.0, 1.5, smoothpar))
    4.01
    """
    return calc_smoothpar_logistic2(metapar)


def calc_smoothpar_min1(metapar):
    """Return the smoothing parameter corresponding to the given meta
    parameter when using |smooth_min1|.

    |smooth_min1| is only an alias for |smooth_logistic2|.

    Calculate the smoothing parameter value corresponding the meta parameter
    value 2.5:

    >>> from hydpy.auxs.smoothtools import calc_smoothpar_min1
    >>> smoothpar = calc_smoothpar_min1(2.5)

    Using this smoothing parameter value, the output of function
    |smooth_min1| is 0.01 below the usual discontinuous minimum
    function result, if the absolute value of the difference
    between the x and the y value is 2.5:

    >>> from hydpy.cythons import smoothutils
    >>> from hydpy import round_
    >>> round_(smoothutils.smooth_min1(-4.0, -1.5, smoothpar))
    -4.01
    """
    return calc_smoothpar_logistic2(metapar)


autodoctools.autodoc_module()

# -*- coding: utf-8 -*-
"""This module implements statistical functionalities frequently used in
hydrological modelling.
"""
# import...
# ...from standard library
from __future__ import division, print_function
# ...from site-packages
import numpy
from hydpy import pandas
from scipy import optimize
from scipy import special
# ...from HydPy
from hydpy.core import autodoctools
from hydpy.core import objecttools
from hydpy.auxs import validtools


def prepare_arrays(sim=None, obs=None, node=None, skip_nan=False):
    """Prepare and return two |numpy| arrays based on the given arguments.

    Note that many functions provided by module |statstools| apply function
    |prepare_arrays| internally (e.g. |nse|).  But you can also apply it
    manually, as shown in the following examples.

    Function |prepare_arrays| can extract time series data from |Node|
    objects.  To set up an example for this, we define a initialization
    time period and prepare a |Node| object:

    >>> from hydpy import pub, Timegrid, Timegrids, Node, round_, nan
    >>> pub.timegrids = Timegrids(Timegrid('01.01.2000',
    ...                                    '07.01.2000',
    ...                                    '1d'))
    >>> node = Node('test')

    Next, we assign values the `simulation` and the `observation` sequences
    (to do so for the `observation` sequence requires a little trick, as
    its values are normally supposed to be read from a file):

    >>> node.prepare_simseries()
    >>> node.sequences.sim.series = 1.0, nan, nan, nan, 2.0, 3.0
    >>> node.sequences.obs.ramflag = True
    >>> node.sequences.obs._setarray([4.0, 5.0, nan, nan, nan, 6.0])

    Now we can pass the node object to function |prepare_arrays| and
    get the (unmodified) time series data:

    >>> from hydpy import prepare_arrays
    >>> arrays = prepare_arrays(node=node)
    >>> round_(arrays[0])
    1.0, nan, nan, nan, 2.0, 3.0
    >>> round_(arrays[1])
    4.0, 5.0, nan, nan, nan, 6.0

    Alternatively, we can pass directly any iterables (e.g. |list| and
    |tuple| objects) containing the `simulated` and `observed` data:

    >>> arrays = prepare_arrays(sim=list(node.sequences.sim.series),
    ...                         obs=tuple(node.sequences.obs.series))
    >>> round_(arrays[0])
    1.0, nan, nan, nan, 2.0, 3.0
    >>> round_(arrays[1])
    4.0, 5.0, nan, nan, nan, 6.0

    The optional `skip_nan` flag allows to skip all values, which are
    no numbers.  Note that only those pairs of `simulated` and `observed`
    values are returned which do not contain any `nan`:

    >>> arrays = prepare_arrays(node=node, skip_nan=True)
    >>> round_(arrays[0])
    1.0, 3.0
    >>> round_(arrays[1])
    4.0, 6.0

    The final examples show the error messages returned in case of
    invalid combinations of input arguments:

    >>> prepare_arrays()
    Traceback (most recent call last):
    ...
    ValueError: Neither a `Node` object is passed to argument `node` nor \
are arrays passed to arguments `sim` and `obs`.

    >>> prepare_arrays(sim=node.sequences.sim.series, node=node)
    Traceback (most recent call last):
    ...
    ValueError: Values are passed to both arguments `sim` and `node`, \
which is not allowed.

    >>> prepare_arrays(obs=node.sequences.obs.series, node=node)
    Traceback (most recent call last):
    ...
    ValueError: Values are passed to both arguments `obs` and `node`, \
which is not allowed.

    >>> prepare_arrays(sim=node.sequences.sim.series)
    Traceback (most recent call last):
    ...
    ValueError: A value is passed to argument `sim` but \
no value is passed to argument `obs`.

    >>> prepare_arrays(obs=node.sequences.obs.series)
    Traceback (most recent call last):
    ...
    ValueError: A value is passed to argument `obs` but \
no value is passed to argument `sim`.
    """
    if node:
        if sim is not None:
            raise ValueError(
                'Values are passed to both arguments `sim` and `node`, '
                'which is not allowed.')
        if obs is not None:
            raise ValueError(
                'Values are passed to both arguments `obs` and `node`, '
                'which is not allowed.')
        sim = node.sequences.sim.series
        obs = node.sequences.obs.series
    elif (sim is not None) and (obs is None):
        raise ValueError(
            'A value is passed to argument `sim` '
            'but no value is passed to argument `obs`.')
    elif (obs is not None) and (sim is None):
        raise ValueError(
            'A value is passed to argument `obs` '
            'but no value is passed to argument `sim`.')
    elif (sim is None) and (obs is None):
        raise ValueError(
            'Neither a `Node` object is passed to argument `node` nor '
            'are arrays passed to arguments `sim` and `obs`.')
    sim = numpy.asarray(sim)
    obs = numpy.asarray(obs)
    if skip_nan:
        idxs = ~numpy.isnan(sim) * ~numpy.isnan(obs)
        sim = sim[idxs]
        obs = obs[idxs]
    return sim, obs


@objecttools.excmessage_decorator(
    'calculate the Nash-Sutcliffe efficiency')
def nse(sim=None, obs=None, node=None, skip_nan=False):
    """Calculate the efficiency criteria after Nash & Sutcliffe.

    If the simulated values predict the observed values as well
    as the average observed value (regarding the the mean square
    error), the NSE value is zero:

    >>> from hydpy import round_
    >>> from hydpy import nse
    >>> nse(sim=[2.0, 2.0, 2.0], obs=[1.0, 2.0, 3.0])
    0.0
    >>> nse(sim=[0.0, 2.0, 4.0], obs=[1.0, 2.0, 3.0])
    0.0

    For worse and better simulated values the NSE is negative
    or positive, respectively:

    >>> nse(sim=[3.0, 2.0, 1.0], obs=[1.0, 2.0, 3.0])
    -3.0
    >>> nse(sim=[1.0, 2.0, 2.0], obs=[1.0, 2.0, 3.0])
    0.5

    The highest possible value is one:

    >>> nse(sim=[1.0, 2.0, 3.0], obs=[1.0, 2.0, 3.0])
    1.0

    See the documentation on function |prepare_arrays| for some
    additional instructions for use of function |nse|.
    """
    sim, obs = prepare_arrays(sim, obs, node, skip_nan)
    return 1.-numpy.sum((sim-obs)**2)/numpy.sum((obs-numpy.mean(obs))**2)


@objecttools.excmessage_decorator(
    'calculate the absolute bias')
def bias_abs(sim=None, obs=None, node=None, skip_nan=False):
    """Calculate the absolute difference between the means of the simulated
    and the observed values.

    >>> from hydpy import round_
    >>> from hydpy import bias_abs
    >>> round_(bias_abs(sim=[2.0, 2.0, 2.0], obs=[1.0, 2.0, 3.0]))
    0.0
    >>> round_(bias_abs(sim=[5.0, 2.0, 2.0], obs=[1.0, 2.0, 3.0]))
    1.0
    >>> round_(bias_abs(sim=[1.0, 1.0, 1.0], obs=[1.0, 2.0, 3.0]))
    -1.0

    See the documentation on function |prepare_arrays| for some
    additional instructions for use of function |bias_abs|.
    """
    sim, obs = prepare_arrays(sim, obs, node, skip_nan)
    return numpy.mean(sim-obs)


@objecttools.excmessage_decorator(
    'calculate the relative bias')
def bias_rel(sim=None, obs=None, node=None, skip_nan=False):
    """Calculate the relative difference between the means of the simulated
    and the observed values.

    >>> from hydpy import round_
    >>> from hydpy import bias_rel
    >>> round_(bias_rel(sim=[2.0, 2.0, 2.0], obs=[1.0, 2.0, 3.0]))
    0.0
    >>> round_(bias_rel(sim=[5.0, 2.0, 2.0], obs=[1.0, 2.0, 3.0]))
    0.5
    >>> round_(bias_rel(sim=[1.0, 1.0, 1.0], obs=[1.0, 2.0, 3.0]))
    -0.5

    See the documentation on function |prepare_arrays| for some
    additional instructions for use of function |bias_rel|.
    """
    sim, obs = prepare_arrays(sim, obs, node, skip_nan)
    return numpy.mean(sim)/numpy.mean(obs)-1.


@objecttools.excmessage_decorator(
    'calculate the standard deviation ratio')
def std_ratio(sim=None, obs=None, node=None, skip_nan=False):
    """Calculate the ratio between the standard deviation of the simulated
    and the observed values.

    >>> from hydpy import round_
    >>> from hydpy import std_ratio
    >>> round_(std_ratio(sim=[1.0, 2.0, 3.0], obs=[1.0, 2.0, 3.0]))
    0.0
    >>> round_(std_ratio(sim=[1.0, 1.0, 1.0], obs=[1.0, 2.0, 3.0]))
    -1.0
    >>> round_(std_ratio(sim=[0.0, 3.0, 6.0], obs=[1.0, 2.0, 3.0]))
    2.0

    See the documentation on function |prepare_arrays| for some
    additional instructions for use of function |std_ratio|.
    """
    sim, obs = prepare_arrays(sim, obs, node, skip_nan)
    return numpy.std(sim)/numpy.std(obs)-1.


@objecttools.excmessage_decorator(
    'calculate the Pearson correlation coefficient')
def corr(sim=None, obs=None, node=None, skip_nan=False):
    """Calculate the product-moment correlation coefficient after Pearson.

    >>> from hydpy import round_
    >>> from hydpy import corr
    >>> round_(corr(sim=[0.5, 1.0, 1.5], obs=[1.0, 2.0, 3.0]))
    1.0
    >>> round_(corr(sim=[4.0, 2.0, 0.0], obs=[1.0, 2.0, 3.0]))
    -1.0
    >>> round_(corr(sim=[1.0, 2.0, 1.0], obs=[1.0, 2.0, 3.0]))
    0.0

    See the documentation on function |prepare_arrays| for some
    additional instructions for use of function |corr|.
    """
    sim, obs = prepare_arrays(sim, obs, node, skip_nan)
    return numpy.corrcoef(sim, obs)[0, 1]


def _pars_sepd(xi, beta):
    # pylint: disable=invalid-name
    gamma1 = special.gamma(3.*(1.+beta)/2.)
    gamma2 = special.gamma((1.+beta)/2.)
    w_beta = gamma1**.5 / (1.+beta) / gamma2**1.5
    c_beta = (gamma1/gamma2)**(1./(1.+beta))
    m_1 = special.gamma(1.+beta) / gamma1**.5 / gamma2**.5
    m_2 = 1.
    mu_xi = m_1*(xi-1./xi)
    sigma_xi = numpy.sqrt((m_2-m_1**2)*(xi**2+1./xi**2)+2*m_1**2-m_2)
    return mu_xi, sigma_xi, w_beta, c_beta


def _pars_h(sigma1, sigma2, sim):
    return sigma1*numpy.mean(sim) + sigma2*sim


@objecttools.excmessage_decorator(
    'calculate the probability densities with the '
    'heteroskedastic skewed exponential power distribution')
def hsepd_pdf(sigma1, sigma2, xi, beta,
              sim=None, obs=None, node=None, skip_nan=False):
    """Calculate the probability densities based on the
    heteroskedastic skewed exponential power distribution.

    For convenience, the required parameters of the probability density
    function as well as the simulated and observed values are stored
    in a dictonary:

    >>> import numpy
    >>> from hydpy import round_
    >>> from hydpy import hsepd_pdf
    >>> general = {'sigma1': 0.2,
    ...            'sigma2': 0.0,
    ...            'xi': 1.0,
    ...            'beta': 0.0,
    ...            'sim': numpy.arange(10.0, 41.0),
    ...            'obs': numpy.full(31, 25.0)}

    The following test function allows the variation of one parameter
    and prints some and plots all of probability density values
    corresponding to different simulated values:

    >>> def test(**kwargs):
    ...     from matplotlib import pyplot
    ...     special = general.copy()
    ...     name, values = list(kwargs.items())[0]
    ...     results = numpy.zeros((len(general['sim']), len(values)+1))
    ...     results[:, 0] = general['sim']
    ...     for jdx, value in enumerate(values):
    ...         special[name] = value
    ...         results[:, jdx+1] = hsepd_pdf(**special)
    ...         pyplot.plot(results[:, 0], results[:, jdx+1],
    ...                     label='%s=%.1f' % (name, value))
    ...     pyplot.legend()
    ...     for idx, result in enumerate(results):
    ...         if not (idx % 5):
    ...             round_(result)

    When varying parameter `beta`, the resulting probabilities correspond
    to the Laplace distribution (1.0), normal distribution (0.0), and the
    uniform distribution (-1.0), respectively.  Note that we use -0.99
    instead of -1.0 for approximating the uniform distribution to prevent
    from running into numerical problems, which are not solved yet:

    >>> test(beta=[1.0, 0.0, -0.99])
    10.0, 0.002032, 0.000886, 0.0
    15.0, 0.008359, 0.010798, 0.0
    20.0, 0.034382, 0.048394, 0.057739
    25.0, 0.141421, 0.079788, 0.057739
    30.0, 0.034382, 0.048394, 0.057739
    35.0, 0.008359, 0.010798, 0.0
    40.0, 0.002032, 0.000886, 0.0

    .. testsetup::

        >>> from matplotlib import pyplot
        >>> pyplot.close()

    When varying parameter `xi`, the resulting density is negatively
    skewed (0.2), symmetric (1.0), and positively skewed (5.0),
    respectively:

    >>> test(xi=[0.2, 1.0, 5.0])
    10.0, 0.0, 0.000886, 0.003175
    15.0, 0.0, 0.010798, 0.012957
    20.0, 0.092845, 0.048394, 0.036341
    25.0, 0.070063, 0.079788, 0.070063
    30.0, 0.036341, 0.048394, 0.092845
    35.0, 0.012957, 0.010798, 0.0
    40.0, 0.003175, 0.000886, 0.0

    .. testsetup::

        >>> from matplotlib import pyplot
        >>> pyplot.close()

    In the above examples, the actual `sigma` (5.0) is calculated by
    multiplying `sigma1` (0.2) with the mean simulated value (25.0),
    internally.  This can be done for modelling homoscedastic errors.
    Instead, `sigma2` is multiplied with the individual simulated values
    to account for heteroscedastic errors.  With increasing values of
    `sigma2`, the resulting densities are modified as follows:

    >>> test(sigma2=[0.0, 0.1, 0.2])
    10.0, 0.000886, 0.002921, 0.005737
    15.0, 0.010798, 0.018795, 0.022831
    20.0, 0.048394, 0.044159, 0.037988
    25.0, 0.079788, 0.053192, 0.039894
    30.0, 0.048394, 0.04102, 0.032708
    35.0, 0.010798, 0.023493, 0.023493
    40.0, 0.000886, 0.011053, 0.015771

    .. testsetup::

        >>> from matplotlib import pyplot
        >>> pyplot.close()
    """
    # pylint: disable=invalid-name
    sim, obs = prepare_arrays(sim, obs, node, skip_nan)
    sigmas = _pars_h(sigma1, sigma2, sim)
    mu_xi, sigma_xi, w_beta, c_beta = _pars_sepd(xi, beta)
    x, mu = obs, sim
    a = (x-mu)/sigmas
    a_xi = numpy.empty(a.shape)
    idxs = mu_xi+sigma_xi*a < 0.
    a_xi[idxs] = numpy.absolute(xi*(mu_xi+sigma_xi*a[idxs]))
    a_xi[~idxs] = numpy.absolute(1./xi*(mu_xi+sigma_xi*a[~idxs]))
    ps = (2.*sigma_xi/(xi+1./xi)*w_beta *
          numpy.exp(-c_beta*a_xi**(2./(1.+beta))))/sigmas
    return ps


def _hsepd_manual(sigma1, sigma2, xi, beta, sim, obs):
    # pylint: disable=invalid-name
    ps = hsepd_pdf(sigma1, sigma2, xi, beta, sim, obs)
    ps[ps < 1e-200] = 1e-200
    return numpy.mean(numpy.log(ps))


@objecttools.excmessage_decorator(
    'calculate an objective value based on method `hsepd_manual`')
def hsepd_manual(sigma1, sigma2, xi, beta,
                 sim=None, obs=None, node=None, skip_nan=False):
    """Calculate the mean of the logarithmised probability densities of the
    'heteroskedastic skewed exponential power distribution.

    The following examples are taken from the documentation of function
    |hsepd_pdf|, which is used by function |hsepd_manual|.  The first
    one deals with a heteroscedastic normal distribution:

    >>> from hydpy import round_
    >>> from hydpy import hsepd_manual
    >>> round_(hsepd_manual(sigma1=0.2, sigma2=0.2,
    ...                     xi=1.0, beta=0.0,
    ...                     sim=numpy.arange(10.0, 41.0),
    ...                     obs=numpy.full(31, 25.0)))
    -3.682842

    The second one is supposed to show to small zero probability density
    values are set to 1e-200 before calculating their logarithm (which
    means that the lowest possible value returned by function
    |hsepd_manual| is approximately -460):

    >>> round_(hsepd_manual(sigma1=0.2, sigma2=0.0,
    ...                     xi=1.0, beta=-0.99,
    ...                     sim=numpy.arange(10.0, 41.0),
    ...                     obs=numpy.full(31, 25.0)))
    -209.539335
    """
    # pylint: disable=invalid-name
    sim, obs = prepare_arrays(sim, obs, node, skip_nan)
    return _hsepd_manual(sigma1, sigma2, xi, beta, sim, obs)


@objecttools.excmessage_decorator(
    'calculate an objective value based on method `hsepd`')
def hsepd(sim=None, obs=None, node=None, skip_nan=False,
          inits=None, return_pars=False, silent=True):
    """Calculate the mean of the logarithmised probability densities of the
    'heteroskedastic skewed exponential power distribution.

    Function |hsepd| serves the same purpose as function |hsepd_manual|,
    but tries to estimate the parameters of the heteroscedastic skewed
    exponential distribution via an optimization algorithm.  This
    is shown by generating a random sample.  1000 simulated values
    are scattered around the observed (true) value of 10.0 with a
    standard deviation of 2.0:

    >>> import numpy
    >>> numpy.random.seed(0)
    >>> sim = numpy.random.normal(10.0, 2.0, 1000)
    >>> obs = numpy.full(1000, 10.0)

    First, as a reference, we calculate the "true" value based on
    function |hsepd_manual| and the correct distribution parameters:

    >>> from hydpy import round_
    >>> from hydpy import hsepd, hsepd_manual
    >>> round_(hsepd_manual(sigma1=0.2, sigma2=0.0,
    ...                     xi=1.0, beta=0.0,
    ...                     sim=sim, obs=obs))
    -2.100093

    When using function |hsepd|, the returned value is even a little
    "better":

    >>> round_(hsepd(sim=sim, obs=obs))
    -2.09983

    This is due to the deviation from the random sample to its
    theoretical distribution.  This is reflected by small differences
    between the estimated values and the theoretical values of
    `sigma1` (0.2), , `sigma2` (0.0), `xi` (1.0), and `beta` (0.0).
    The estimated values are returned in the mentioned order through
    enabling the `return_pars` option:

    >>> value, pars = hsepd(sim=sim, obs=obs, return_pars=True)
    >>> round_(pars, decimals=5)
    0.19966, 0.0, 0.96836, 0.0188

    There is no guarantee that the optimization numerical optimization
    algorithm underlying function |hsepd| will always find the parameters
    resulting in the largest value returned by function |hsepd_manual|.
    You can increase its robustness (and decrease computation time) by
    supplying good initial parameter values:

    >>> value, pars = hsepd(sim=sim, obs=obs, return_pars=True,
    ...                     inits=(0.2, 0.0, 1.0, 0.0))
    >>> round_(pars, decimals=5)
    0.19966, 0.0, 0.96836, 0.0188

    However, the following example shows a case when this strategie
    results in worse results:

    >>> value, pars = hsepd(sim=sim, obs=obs, return_pars=True,
    ...                     inits=(0.0, 0.2, 1.0, 0.0))
    >>> round_(value)
    -2.174492
    >>> round_(pars)
    0.0, 0.213179, 1.705485, 0.505112
    """
    # pylint: disable=invalid-name

    def transform(pars):
        """Transform the actual optimization problem into a function to
        be minimized and apply parameter constraints."""
        sigma1, sigma2, xi, beta = constrain(*pars)
        return -_hsepd_manual(sigma1, sigma2, xi, beta, sim, obs)

    def constrain(sigma1, sigma2, xi, beta):
        """Apply constrains on the given parameter values."""
        sigma1 = numpy.clip(sigma1, 0.0, None)
        sigma2 = numpy.clip(sigma2, 0.0, None)
        xi = numpy.clip(xi, 0.1, 10.0)
        beta = numpy.clip(beta, -0.99, 5.0)
        return sigma1, sigma2, xi, beta

    sim, obs = prepare_arrays(sim, obs, node, skip_nan)
    if not inits:
        inits = [0.1, 0.2, 3.0, 1.0]
    values = optimize.fmin(transform, inits,
                           ftol=1e-12, xtol=1e-12,
                           disp=not silent)
    values = constrain(*values)   # pylint: disable=too-many-function-args
    result = _hsepd_manual(*values, sim=sim, obs=obs)
    if return_pars:
        return result, values
    return result


@objecttools.excmessage_decorator(
    'calculate the weighted mean time')
def calc_mean_time(timepoints, weights):
    """Return the weighted mean of the given timepoints.

    With equal given weights, the result is simply the mean of the given
    time points:

    >>> from hydpy import calc_mean_time
    >>> calc_mean_time(timepoints=[3., 7.],
    ...                weights=[2., 2.])
    5.0

    With different weights, the resulting mean time is shifted to the larger
    ones:

    >>> calc_mean_time(timepoints=[3., 7.],
    ...                weights=[1., 3.])
    6.0

    Or, in the most extreme case:

    >>> calc_mean_time(timepoints=[3., 7.],
    ...                weights=[0., 4.])
    7.0

    There will be some checks for input plausibility perfomed, e.g.:

    >>> calc_mean_time(timepoints=[3., 7.],
    ...                weights=[-2., 2.])
    Traceback (most recent call last):
    ...
    ValueError: While trying to calculate the weighted mean time, \
the following error occured: For the following objects, at least \
one value is negative: weights.
    """
    timepoints = numpy.array(timepoints)
    weights = numpy.array(weights)
    validtools.test_equal_shape(timepoints=timepoints, weights=weights)
    validtools.test_non_negative(weights=weights)
    return numpy.dot(timepoints, weights)/numpy.sum(weights)


@objecttools.excmessage_decorator(
    'calculate the weighted time deviation from mean time')
def calc_mean_time_deviation(timepoints, weights, mean_time=None):
    """Return the weighted deviation of the given timepoints from their mean
    time.

    With equal given weights, the is simply the standard deviation of the
    given time points:

    >>> from hydpy import calc_mean_time_deviation
    >>> calc_mean_time_deviation(timepoints=[3., 7.],
    ...                          weights=[2., 2.])
    2.0

    One can pass a precalculated or alternate mean time:

    >>> from hydpy import round_
    >>> round_(calc_mean_time_deviation(timepoints=[3., 7.],
    ...                                 weights=[2., 2.],
    ...                                 mean_time=4.))
    2.236068

    >>> round_(calc_mean_time_deviation(timepoints=[3., 7.],
    ...                                 weights=[1., 3.]))
    1.732051

    Or, in the most extreme case:

    >>> calc_mean_time_deviation(timepoints=[3., 7.],
    ...                          weights=[0., 4.])
    0.0

    There will be some checks for input plausibility perfomed, e.g.:

    >>> calc_mean_time_deviation(timepoints=[3., 7.],
    ...                          weights=[-2., 2.])
    Traceback (most recent call last):
    ...
    ValueError: While trying to calculate the weighted time deviation \
from mean time, the following error occured: For the following objects, \
at least one value is negative: weights.
    """
    timepoints = numpy.array(timepoints)
    weights = numpy.array(weights)
    validtools.test_equal_shape(timepoints=timepoints, weights=weights)
    validtools.test_non_negative(weights=weights)
    if mean_time is None:
        mean_time = calc_mean_time(timepoints, weights)
    return (numpy.sqrt(numpy.dot(weights, (timepoints-mean_time)**2) /
                       numpy.sum(weights)))


@objecttools.excmessage_decorator(
    'evaluate the simulation results of some node objects')
def evaluationtable(nodes, criteria, nodenames=None,
                    critnames=None, skip_nan=False):
    """Return a table containing the results of the given evaluation
    criteria for the given |Node| objects.

    First, we define two nodes with different simulation and observation
    data (see function |prepare_arrays| for some explanations):

    >>> from hydpy import pub, Timegrid, Timegrids, Node, round_, nan
    >>> pub.timegrids = Timegrids(Timegrid('01.01.2000',
    ...                                    '04.01.2000',
    ...                                    '1d'))
    >>> nodes = Node('test1'), Node('test2')
    >>> for node in nodes:
    ...     node.prepare_simseries()
    ...     node.sequences.sim.series = 1.0, 2.0, 3.0
    ...     node.sequences.obs.ramflag = True
    ...     node.sequences.obs._setarray([4.0, 5.0, 6.0])
    >>> nodes[0].sequences.sim.series = 1.0, 2.0, 3.0
    >>> nodes[0].sequences.obs._setarray([4.0, 5.0, 6.0])
    >>> nodes[1].sequences.sim.series = 1.0, 2.0, 3.0
    >>> nodes[1].sequences.obs._setarray([3.0, nan, 1.0])

    Selecting functions |corr| and |bias_abs| as evaluation criteria,
    function |evaluationtable| returns the following table (which is
    actually a pandas data frame):

    >>> from hydpy import evaluationtable, corr, bias_abs
    >>> evaluationtable(nodes, (corr, bias_abs))
           corr  bias_abs
    test1   1.0      -3.0
    test2   NaN       NaN

    One can pass alternative names for both the node objects and the
    criteria functions.  Also, `nan` values can be skipped:

    >>> evaluationtable(nodes, (corr, bias_abs),
    ...                 nodenames=('first node', 'second node'),
    ...                 critnames=('corrcoef', 'bias'),
    ...                 skip_nan=True)
                 corrcoef  bias
    first node        1.0  -3.0
    second node      -1.0   0.0

    The number of assigned node objects and criteria functions must
    match the number of givern alternative names:

    >>> evaluationtable(nodes, (corr, bias_abs),
    ...                 nodenames=('first node',))
    Traceback (most recent call last):
    ...
    ValueError: While trying to evaluate the simulation results of some \
node objects, the following error occured: 2 node objects are given \
which does not match with number of given alternative names beeing 1.

    >>> evaluationtable(nodes, (corr, bias_abs),
    ...                 critnames=('corrcoef',))
    Traceback (most recent call last):
    ...
    ValueError: While trying to evaluate the simulation results of some \
node objects, the following error occured: 2 criteria functions are given \
which does not match with number of given alternative names beeing 1.
    """
    if nodenames:
        if len(nodes) != len(nodenames):
            raise ValueError(
                '%d node objects are given which does not match with '
                'number of given alternative names beeing %s.'
                % (len(nodes), len(nodenames)))
    else:
        nodenames = [node.name for node in nodes]
    if critnames:
        if len(criteria) != len(critnames):
            raise ValueError(
                '%d criteria functions are given which does not match '
                'with number of given alternative names beeing %s.'
                % (len(criteria), len(critnames)))
    else:
        critnames = [crit.__name__ for crit in criteria]
    data = numpy.empty((len(nodes), len(criteria)), dtype=float)
    for idx, node in enumerate(nodes):
        sim, obs = prepare_arrays(None, None, node, skip_nan)
        for jdx, criterion in enumerate(criteria):
            data[idx, jdx] = criterion(sim, obs)
    table = pandas.DataFrame(
        data=data, index=nodenames, columns=critnames)
    return table


autodoctools.autodoc_module()

# import...
"""This module implements features for the validation of (numerical) input
data.
"""
# ...from standard library
from __future__ import division, print_function
# ...from site-packages
import numpy
# ...from HydPy
from hydpy.core import autodoctools


def test_equal_shape(**kwargs):
    """Raise a ValueError if the shapes of the objects given as keywords
    are not equal.

    If all shapes are equal, nothing happens:

    >>> from hydpy.auxs.validtools import test_equal_shape
    >>> test_equal_shape(arr1=numpy.array([1., 2.]),
    ...                  arr2=numpy.array([3., 4.]),
    ...                  arr3=numpy.array([5., 6.]))

    If at least one shape differs, the following error is raised:

    >>> test_equal_shape(arr1=numpy.array([1., 2.]),
    ...                  arr2=numpy.array([3.]),
    ...                  arr3=numpy.array([5., 6.]))
    Traceback (most recent call last):
    ...
    ValueError: The shapes of the following objects are not equal: arr1 (2,), arr2 (1,), arr3 (2,).

    For flexibility in the functions application, it is allowed to pass only
    one array or no arrays at all:

    >>> test_equal_shape(arr1=numpy.array([1., 2.]))
    >>> test_equal_shape()
    """
    names = list(kwargs.keys())
    shapes = numpy.array([numpy.array(array).shape
                          for array in kwargs.values()])
    if any(shapes[:-1] != shapes[1:]):
        raise ValueError(
            'The shapes of the following objects are not equal: %s.'
            % ', '.join('%s %s' % (name, tuple(shape)) for (name, shape)
                        in sorted(zip(names, shapes))))


def test_non_negative(**kwargs):
    """Raise a ValueError if at least one value of the objects given as
    keywords is negative.

    If all values are non negative, nothing happens:

    >>> from hydpy.auxs.validtools import test_non_negative
    >>> test_non_negative(arr1=numpy.array([1., 2.]),
    ...                   arr2=numpy.array([3., 4.]),
    ...                   arr3=numpy.array([5., 6.]))

    If at least one value is negative, the following error is raised:

    >>> test_non_negative(arr1=numpy.array([1., 2.]),
    ...                   arr2=numpy.array([-3., 4.]),
    ...                   arr3=numpy.array([5., 6.]))
    Traceback (most recent call last):
    ...
    ValueError: For the following objects, at least one value is negative: arr2.

    For flexibility in the functions application, it is allowed to pass
    no arrays at all:

    >>> test_non_negative()
    """
    names = list(kwargs.keys())
    negs = [numpy.nanmin(array) < 0. for array in kwargs.values()]
    if any(negs):
        raise ValueError(
            'For the following objects, at least one value is negative: %s.'
            % ', '.join(name for name, neg in sorted(zip(names, negs)) if neg))

autodoctools.autodoc_module()

# -*- coding: utf-8 -*-
"""This module provides some abstract base classes.

There are some type checks within the HydPy framework relying on the
build in  function |isinstance|.  In order to keep HydPy "pythonic",
the following abstract base classes are defined.  All calls to |isinstance|
should rely these abstract base classes instead of the respective
original classes.  This helps to build e.g. a new parameter class when
one wants to avoid to inherit from |Parameter|.

At the moment, the provided classes do not provide things like abstract
methods (should be added later).  Just use them to register new classes
that are not actual subclasses of the respective HydPy classes, but
should be handled as if they were.  See class |anntools.ANN| as an example.
"""
# import...
# ...from standard library
from __future__ import division, print_function
import abc
# ...from HydPy
from hydpy import pub
from hydpy.core import autodoctools


if pub.pyversion > 2:
    _ABC = abc.ABC
else:
    class _ABC(object):
        __metaclass__ = abc.ABCMeta


class DocABC(_ABC):
    """ABC base class automatically documenting is registered subclasses."""

    _registry_empty = True

    @classmethod
    def register(cls, subclass):
        """Add information to the documentation of the given abstract base
        class and register the subclass afterwards.

        Subclass the new abstract base class `NewABC` and define some new
        concrete classes (`New1`, `New2`, `New3`) which do not inherit
        from `NewABC`:

        >>> from hydpy.core.abctools import DocABC
        >>> class NewABC(DocABC):
        ...    "A new base class."
        >>> class New1(object):
        ...     "First new class"
        >>> class New2(object):
        ...     "Second new class"
        >>> class New3(object):
        ...     "Third new class"

        The docstring `NewABC` is still the same as defined above:

        >>> print(NewABC.__doc__)
        A new base class.

        Now we register the concrete classes `New1` and `New2`:

        >>> NewABC.register(New2)
        >>> NewABC.register(New1)
        >>> NewABC.register(New2)

        Now the docstring of `NewABC` includes the information about
        the concrete classes already registered:

        >>> print(NewABC.__doc__)
        A new base class.
        <BLANKLINE>
        At the moment, the following classes are registered:
             * :class:`~hydpy.core.abctools.New2`
             * :class:`~hydpy.core.abctools.New1`

        Note that the docstring order is the registration order.
        Also note that the "accidental reregistration" of class
        `New2` does not modify the docstring.

        Now the concrete classes `New1` and `New2` are handled as
        if they were actual subclasses of `NewABC`, but class `New3`
        -- which had not been registered -- is not:

        >>> issubclass(New1, NewABC)
        True
        >>> isinstance(New1(), NewABC)
        True
        >>> issubclass(New2, NewABC)
        True
        >>> isinstance(New2(), NewABC)
        True
        >>> issubclass(New3, NewABC)
        False
        >>> isinstance(New3(), NewABC)
        False
        """
        if cls._registry_empty:
            cls._registry_empty = False
            cls.__doc__ += \
                '\n\nAt the moment, the following classes are registered:'
        if not issubclass(subclass, cls):
            cls.__doc__ += ('\n     * :class:`~%s`'
                            % str(subclass).split("'")[1])
            abc.ABCMeta.register(cls, subclass)


class IterableNonStringABC(_ABC):
    """Abstract base class for checking if an object is iterable but not a
    string."""

    @classmethod
    def __subclasshook__(cls, C):
        if cls is IterableNonStringABC:
            return (hasattr(C, '__iter__') and
                    not (isinstance(C, str) or
                         issubclass(C, str)))
        return NotImplemented


class ElementABC(DocABC):
    """Abstract base class for registering custom |Element| classes."""
    pass


class NodeABC(DocABC):
    """Abstract base class for registering custom |Node| classes."""
    pass


class VariableABC(DocABC):
    """Abstract base class for registering custom |Variable| classes.

    Usually, new classes should either be registered as a parameter
    or a sequence.  Afterwards, they are automatically handled as
    |Variable| subclasses:

    >>> from hydpy.core.abctools import VariableABC, ParameterABC
    >>> class New(object):
    ...     pass
    >>> issubclass(New, VariableABC)
    False
    >>> ParameterABC.register(New)
    >>> issubclass(New, VariableABC)
    True
    """


class ParameterABC(VariableABC):
    """Abstract base class for registering custom |Parameter| classes."""


class ANNABC(DocABC):
    """Abstract base class for registering custom |anntools.ANN| classes."""


class SeasonalANNABC(DocABC):
    """Abstract base class for registering custom |anntools.SeasonalANN|
    classes."""


class IOSequencesABC(DocABC):
    """Abstract base class for registering custom |IOSequences| classes."""


class InputSequencesABC(DocABC):
    """Abstract base class for registering custom |InputSequences| classes."""


class OutputSequencesABC(DocABC):
    """Abstract base class for registering custom "OutputSequences" classes
    like |FluxSequences|."""


class SequenceABC(VariableABC):
    """Abstract base class for registering custom |Sequence| classes."""
    pass


class InputSequenceABC(SequenceABC):
    """Abstract base class for registering custom |InputSequence| classes."""
    pass


class FluxSequenceABC(SequenceABC):
    """Abstract base class for registering custom |FluxSequence| classes."""
    pass


class ConditionSequenceABC(SequenceABC):
    """Abstract base class for registering custom |ConditionSequence| classes.
    """
    pass


class StateSequenceABC(ConditionSequenceABC):
    """Abstract base class for registering custom |StateSequence| classes."""
    pass


class LogSequenceABC(ConditionSequenceABC):
    """Abstract base class for registering custom |LogSequence| classes."""
    pass


class AideSequenceABC(SequenceABC):
    """Abstract base class for registering custom |AideSequence| classes."""
    pass


class LinkSequenceABC(SequenceABC):
    """Abstract base class for registering custom |LinkSequence| classes."""
    pass


class NodeSequenceABC(SequenceABC):
    """Abstract base class for registering custom |NodeSequence| classes."""
    pass


class DateABC(DocABC):
    """Abstract base class for registering custom |Date| classes."""
    pass


class PeriodABC(DocABC):
    """Abstract base class for registering custom |Period| classes."""
    pass


class TimegridABC(DocABC):
    """Abstract base class for registering custom |Timegrid| classes."""
    pass


class TimegridsABC(DocABC):
    """Abstract base class for registering custom |Timegrids| classes."""
    pass


class TOYABC(DocABC):
    """Abstract base class for registering custom |TOY| classes."""
    pass


class ModelABC(DocABC):
    """Abstract base class for registering custom |Model| classes."""
    pass


autodoctools.autodoc_module()

# -*- coding: utf-8 -*-
"""This module implements tools for increasing the level of automation and
standardisation of the online documentation generated with Sphinx.
"""
# import...
# ...from the Python standard library
from __future__ import division, print_function
import collections
import copy
import datetime
import doctest
import importlib
import inspect
import io
import os
import pkgutil
import types
import unittest
# ...from site-packages
import numpy
import wrapt
# ...from HydPy
import hydpy
from hydpy import builtins
from hydpy import pub
from hydpy import config
from hydpy import auxs
from hydpy import core
from hydpy import cythons
from hydpy import models
from hydpy.cythons.autogen import annutils
from hydpy.cythons.autogen import pointerutils
from hydpy.cythons.autogen import smoothutils


def description(self):
    """Returns the first "paragraph" of the docstring of the given object.

    Note that ugly things like multiple whitespaces and newline characters
    are removed:

    >>> from hydpy.core import autodoctools, objecttools
    >>> autodoctools.description(objecttools.augment_excmessage)
    'Augment an exception message with additional information while keeping \
the original traceback.'

    In case the given object does not define a docstring, the following
    is returned:
    >>> autodoctools.description(type('Test', (), {}))
    'no description available'
    """
    if self.__doc__ in (None, ''):
        return 'no description available'
    return ' '.join(self.__doc__.split('\n\n')[0].split())


_PAR_SPEC2CAPT = collections.OrderedDict((('parameters', 'Parameter tools'),
                                          ('constants', 'Constants'),
                                          ('control', 'Control parameters'),
                                          ('derived', 'Derived parameters')))

_SEQ_SPEC2CAPT = collections.OrderedDict((('sequences', 'Sequence tools'),
                                          ('inputs', 'Input sequences'),
                                          ('fluxes', 'Flux sequences'),
                                          ('states', 'State sequences'),
                                          ('logs', 'Log sequences'),
                                          ('inlets', 'Inlet sequences'),
                                          ('outlets', 'Outlet sequences'),
                                          ('receivers', 'Receiver sequences'),
                                          ('senders', 'Sender sequences'),
                                          ('aides', 'Aide sequences')))

_all_spec2capt = _PAR_SPEC2CAPT.copy()   # pylint: disable=invalid-name
_all_spec2capt.update(_SEQ_SPEC2CAPT)


@wrapt.decorator
def make_autodoc_optional(wrapped, instance, args, kwargs):
    """Decorate function related to automatic documentation refinement,
    so that they will be applied only when requested (when `use_autodoc`
    of module |config| is `True`) or when possible (when `HydPy` is not
    frozen/bundled)."""
    if config.use_autodoc and not pub._is_hydpy_bundled:
        return wrapped(*args, **kwargs)
    return None


def _add_title(title, marker):
    """Return a title for a basemodels docstring."""
    return ['', title, marker*len(title)]


def _add_lines(specification, module):
    """Return autodoc commands for a basemodels docstring.

    Note that `collection classes` (e.g. `Model`, `ControlParameters`,
    `InputSequences` are placed on top of the respective section and the
    `contained classes` (e.g. model methods, `ControlParameter` instances,
    `InputSequence` instances at the bottom.  This differs from the order
    of their definition in the respective modules, but results in a better
    documentation structure.
    """
    caption = _all_spec2capt.get(specification, 'dummy')
    if caption.split()[-1] in ('parameters', 'sequences'):
        exists_collectionclass = True
        name_collectionclass = caption.title().replace(' ', '')
    else:
        exists_collectionclass = False
    lines = []
    if specification == 'model':
        lines += ['',
                  '.. autoclass:: ' + module.__name__ + '.Model',
                  '    :members:',
                  '    :show-inheritance:']
    elif exists_collectionclass:
        lines += ['',
                  '.. autoclass:: %s.%s' % (module.__name__,
                                            name_collectionclass),
                  '    :members:',
                  '    :show-inheritance:']
    lines += ['',
              '.. automodule:: ' + module.__name__,
              '    :members:',
              '    :show-inheritance:']
    if specification == 'model':
        lines += ['    :exclude-members: Model']
    elif exists_collectionclass:
        lines += ['    :exclude-members: ' + name_collectionclass]
    return lines


@make_autodoc_optional
def autodoc_basemodel():
    """Add an exhaustive docstring to the `__init__` module of a basemodel.

    One just has to write `autodoc_basemodel()` at the bottom of an `__init__`
    module of a basemodel, and all model, parameter and sequence information
    are appended to the modules docstring.  The resulting docstring is suitable
    automatic documentation generation via `Sphinx` and `autodoc`.  Hence
    it helps in constructing HydPy's online documentation and supports the
    embeded help feature of `Spyder` (to see the result, import the package
    of an arbitrary basemodel, e.g. `from hydpy.models import lland` and
    press `cntr+i` with the cursor placed on `lland` written in the IPython
    console afterwards).

    Note that the resulting documentation will be complete only when the
    modules of the basemodel are named in the standard way, e.g. `lland_model`,
    `lland_control`, `lland_inputs`.
    """
    namespace = inspect.currentframe().f_back.f_back.f_locals
    doc = namespace.get('__doc__')
    if doc is None:
        doc = ''
    basemodulename = namespace['__name__'].split('.')[-1]
    modules = {key: value for key, value in namespace.items()
               if (isinstance(value, types.ModuleType) and
                   key.startswith(basemodulename+'_'))}
    substituter = Substituter(hydpy.substituter)
    lines = []
    specification = 'model'
    modulename = basemodulename+'_'+specification
    if modulename in modules:
        module = modules[modulename]
        lines += _add_title('Model features', '-')
        lines += _add_lines(specification, module)
        substituter.add_module(module)
    for (spec2capt, title) in zip((_PAR_SPEC2CAPT, _SEQ_SPEC2CAPT),
                                  ('Parameter features', 'Sequence features')):
        new_lines = _add_title(title, '-')
        found_module = False
        for (specification, caption) in spec2capt.items():
            modulename = basemodulename+'_'+specification
            module = modules.get(modulename)
            if module:
                found_module = True
                new_lines += _add_title(caption, '.')
                new_lines += _add_lines(specification, module)
                substituter.add_module(module)
        if found_module:
            lines += new_lines
    doc += '\n'.join(lines)
    namespace['__doc__'] = doc
    basemodule = importlib.import_module(namespace['__name__'])
    substituter.add_module(basemodule)
    substituter.update_masters()
    namespace['substituter'] = substituter


@make_autodoc_optional
def autodoc_applicationmodel():
    """Improves the docstrings of application models when called
    at the bottom of the respective module.

    |autodoc_applicationmodel| requires, similar to
    |autodoc_basemodel|, that both the application model and its
    base model are defined in the conventional way.
    """
    namespace = inspect.currentframe().f_back.f_back.f_locals
    doc = namespace.get('__doc__')
    if doc is None:
        doc = ''
    name_applicationmodel = namespace['__name__']
    module_applicationmodel = importlib.import_module(name_applicationmodel)
    name_basemodel = name_applicationmodel.split('_')[0]
    module_basemodel = importlib.import_module(name_basemodel)
    substituter = Substituter(module_basemodel.substituter)
    substituter.add_module(module_applicationmodel)
    substituter.update_masters()
    namespace['substituter'] = substituter


class Substituter(object):
    """Implements a HydPy specific docstring substitution mechanism."""

    def __init__(self, master=None):
        self.master = master
        if master:
            self._short2long = copy.deepcopy(master._short2long)
            self._medium2long = copy.deepcopy(master._medium2long)
            self._blacklist = copy.deepcopy(master._blacklist)
        else:
            self._short2long = {}
            self._medium2long = {}
            self._blacklist = set()

    @staticmethod
    def consider_member(name_member, member, module):
        """Return |True| if the given member should be add to the
        substitutions. If not return |False|.

        Some examples based on the site-package |numpy|:

        >>> from hydpy.core.autodoctools import Substituter
        >>> import numpy

        A constant like |nan| should be added:

        >>> Substituter.consider_member(
        ...     'nan', numpy.nan, numpy)
        True

        Members with a prefixed underscore should not be added:

        >>> Substituter.consider_member(
        ...     '_NoValue', numpy._NoValue, numpy)
        False

        Members that are actually imported modules should not be added:

        >>> Substituter.consider_member(
        ...     'warnings', numpy.warnings, numpy)
        False

        Members that are actually defined in other modules should
        not be added:

        >>> numpy.Substituter = Substituter
        >>> Substituter.consider_member(
        ...     'Substituter', numpy.Substituter, numpy)
        False
        >>> del numpy.Substituter

        Members that are defined in submodules of a given package
        (either from the standard library or from site-packages)
        should be added...

        >>> Substituter.consider_member(
        ...     'clip', numpy.clip, numpy)
        True

        ...but not members defined in :ref:`HydPy` submodules:

        >>> import hydpy
        >>> Substituter.consider_member(
        ...     'Node', hydpy.Node, hydpy)
        False
        """
        if name_member.startswith('_'):
            return False
        if inspect.ismodule(member):
            return False
        real_module = getattr(member, '__module__', None)
        if not real_module:
            return True
        if real_module != module.__name__:
            if 'hydpy' in real_module:
                return False
            elif module.__name__ not in real_module:
                return False
        return True

    @staticmethod
    def get_role(member, cython=False):
        """Return the reStructuredText role `func`, `class`, or `const`
        best describing the given member.

        Some examples based on the site-package |numpy|.  |numpy.clip|
        is a function:

        >>> from hydpy.core.autodoctools import Substituter
        >>> import numpy
        >>> Substituter.get_role(numpy.clip)
        'func'

        |numpy.ndarray| is a class:

        >>> Substituter.get_role(numpy.ndarray)
        'class'

        |numpy.ndarray.clip| is a method, for which also the `function`
        role is returned:

        >>> Substituter.get_role(numpy.ndarray.clip)
        'func'

        For everything else the `constant` role is returned:

        >>> Substituter.get_role(numpy.nan)
        'const'

        When analysing cython extension modules, set the option `cython`
        flag to |True|.  |Double| is correctly identified as a class:

        >>> from hydpy.cythons import pointerutils
        >>> Substituter.get_role(pointerutils.Double, cython=True)
        'class'

        Only with the `cython` flag beeing |True|, for everything else
        the `function` text role is returned (doesn't make sense here,
        but the |numpy| module is not something defined in module
        |pointerutils| anyway):

        >>> Substituter.get_role(pointerutils.numpy, cython=True)
        'func'
        """
        if inspect.isroutine(member) or isinstance(member, numpy.ufunc):
            return 'func'
        elif inspect.isclass(member):
            return 'class'
        elif cython:
            return 'func'
        return 'const'

    def add_substitution(self, short, medium, long, module):
        """Add the given substitutions both as a `short2long` and a
        `medium2long` mapping.

        Assume `variable1` is defined in the hydpy module `module1` and the
        short and medium descriptions are `var1` and `mod1.var1`:

        >>> import types
        >>> module1 = types.ModuleType('hydpy.module1')
        >>> from hydpy.core.autodoctools import Substituter
        >>> substituter = Substituter()
        >>> substituter.add_substitution(
        ...     'var1', 'mod1.var1', 'module1.variable1', module1)
        >>> print(substituter.get_commands())
        .. var1 replace:: module1.variable1
        .. mod1.var1 replace:: module1.variable1

        Adding `variable2` of `module2` has no effect on the predefined
        substitutions:

        >>> module2 = types.ModuleType('hydpy.module2')
        >>> substituter.add_substitution(
        ...     'var2', 'mod2.var2', 'module2.variable2', module2)
        >>> print(substituter.get_commands())
        .. var1 replace:: module1.variable1
        .. var2 replace:: module2.variable2
        .. mod1.var1 replace:: module1.variable1
        .. mod2.var2 replace:: module2.variable2

        But when adding `variable1` of `module2`, the `short2long` mapping
        of `variable1` would become inconclusive, which is why the new
        one (related to `module2`) is not stored and the old one (related
        to `module1`) is removed:

        >>> substituter.add_substitution(
        ...     'var1', 'mod2.var1', 'module2.variable1', module2)
        >>> print(substituter.get_commands())
        .. var2 replace:: module2.variable2
        .. mod1.var1 replace:: module1.variable1
        .. mod2.var1 replace:: module2.variable1
        .. mod2.var2 replace:: module2.variable2

        Adding `variable2` of `module2` accidentally again, does not
        result in any undesired side-effects:

        >>> substituter.add_substitution(
        ...     'var2', 'mod2.var2', 'module2.variable2', module2)
        >>> print(substituter.get_commands())
        .. var2 replace:: module2.variable2
        .. mod1.var1 replace:: module1.variable1
        .. mod2.var1 replace:: module2.variable1
        .. mod2.var2 replace:: module2.variable2

        In order to reduce the risk of name conflicts, only the
        `medium2long` mapping is supported for modules not part of the
        :ref:`HydPy` package:

        >>> module3 = types.ModuleType('module3')
        >>> substituter.add_substitution(
        ...     'var3', 'mod3.var3', 'module3.variable3', module3)
        >>> print(substituter.get_commands())
        .. var2 replace:: module2.variable2
        .. mod1.var1 replace:: module1.variable1
        .. mod2.var1 replace:: module2.variable1
        .. mod2.var2 replace:: module2.variable2
        .. mod3.var3 replace:: module3.variable3

        The only exception to this rule is |builtins|, for which only
        the `short2long` mapping is supported (note also, that the
        module name `builtins` is removed from string `long`):

        >>> from hydpy import builtins
        >>> substituter.add_substitution(
        ...     'str', 'blt.str', ':func:`~builtins.str`', builtins)
        >>> print(substituter.get_commands())
        .. str replace:: :func:`str`
        .. var2 replace:: module2.variable2
        .. mod1.var1 replace:: module1.variable1
        .. mod2.var1 replace:: module2.variable1
        .. mod2.var2 replace:: module2.variable2
        .. mod3.var3 replace:: module3.variable3
        """
        name = module.__name__
        if 'builtin' in name:
            self._short2long[short] = long.split('~')[0] + long.split('.')[-1]
        else:
            if ('hydpy' in name) and (short not in self._blacklist):
                if short in self._short2long:
                    if self._short2long[short] != long:
                        self._blacklist.add(short)
                        del self._short2long[short]
                else:
                    self._short2long[short] = long
            self._medium2long[medium] = long

    def add_module(self, module, cython=False):
        """Add the given module, its members, and their submembers.

        The first examples are based on the site-package |numpy|: which
        is passed to method |Substituter.add_module|:

        >>> from hydpy.core.autodoctools import Substituter
        >>> substituter = Substituter()
        >>> import numpy
        >>> substituter.add_module(numpy)

        Firstly, the module itself is added:

        >>> substituter.find('|numpy|')
        |numpy| :mod:`~numpy`

        Secondly, constants like |numpy.nan| are added:

        >>> substituter.find('|numpy.nan|')
        |numpy.nan| :const:`~numpy.nan`

        Thirdly, functions like |numpy.clip| are added:

        >>> substituter.find('|numpy.clip|')
        |numpy.clip| :func:`~numpy.clip`

        Fourthly, clases line |numpy.ndarray| are added:

        >>> substituter.find('|numpy.ndarray|')
        |numpy.ndarray| :class:`~numpy.ndarray`

        When adding Cython modules, the `cython` flag should be set |True|:

        >>> from hydpy.cythons import pointerutils
        >>> substituter.add_module(pointerutils, cython=True)
        >>> substituter.find('set_pointer')
        |PPDouble.set_pointer| \
:func:`~hydpy.cythons.autogen.pointerutils.PPDouble.set_pointer`
        |pointerutils.PPDouble.set_pointer| \
:func:`~hydpy.cythons.autogen.pointerutils.PPDouble.set_pointer`
        """
        name_module = module.__name__.split('.')[-1]
        short = ('|%s|'
                 % name_module)
        long = (':mod:`~%s`'
                % module.__name__)
        self._short2long[short] = long
        for (name_member, member) in vars(module).items():
            if self.consider_member(
                    name_member, member, module):
                role = self.get_role(member, cython)
                short = ('|%s|'
                         % name_member)
                medium = ('|%s.%s|'
                          % (name_module,
                             name_member))
                long = (':%s:`~%s.%s`'
                        % (role,
                           module.__name__,
                           name_member))
                self.add_substitution(short, medium, long, module)
                if inspect.isclass(member):
                    for name_submember, submember in vars(member).items():
                        if self.consider_member(
                                name_submember, submember, module):
                            role = self.get_role(submember, cython)
                            short = ('|%s.%s|'
                                     % (name_member,
                                        name_submember))
                            medium = ('|%s.%s.%s|'
                                      % (name_module,
                                         name_member,
                                         name_submember))
                            long = (':%s:`~%s.%s.%s`'
                                    % (role,
                                       module.__name__,
                                       name_member,
                                       name_submember))
                            self.add_substitution(short, medium, long, module)

    def add_modules(self, package):
        """Add the modules of the given package without their members."""
        for name in os.listdir(package.__path__[0]):
            if name.startswith('_'):
                continue
            name = name.split('.')[0]
            short = '|%s|' % name
            long = ':mod:`~%s.%s`' % (package.__package__, name)
            self._short2long[short] = long

    def update_masters(self):
        """Update all `master` |Substituter| objects.

        If a |Substituter| objects is passed to the constructor of another
        |Substituter| object, it becomes its `master`:

        >>> from hydpy.core.autodoctools import Substituter
        >>> sub1 = Substituter()
        >>> from hydpy.core import devicetools
        >>> sub1.add_module(devicetools)
        >>> sub2 = Substituter(sub1)
        >>> sub3 = Substituter(sub2)
        >>> sub3.master.master is sub1
        True

        During initialization, all mappings handled by the master object
        are passed to its slave:

        >>> sub3.find('Node|')
        |Node| :class:`~hydpy.core.devicetools.Node`
        |devicetools.Node| :class:`~hydpy.core.devicetools.Node`


        Updating a slave, does not affect its master directly:

        >>> from hydpy.core import hydpytools
        >>> sub3.add_module(hydpytools)
        >>> sub3.find('HydPy|')
        |HydPy| :class:`~hydpy.core.hydpytools.HydPy`
        |hydpytools.HydPy| :class:`~hydpy.core.hydpytools.HydPy`
        >>> sub2.find('HydPy|')

        Through calling |Substituter.update_masters|, the `medium2long`
        mappings are passed the slave to its master:

        >>> sub3.update_masters()
        >>> sub2.find('HydPy|')
        |hydpytools.HydPy| :class:`~hydpy.core.hydpytools.HydPy`

        Then each master object updates its own master object also:

        >>> sub1.find('HydPy|')
        |hydpytools.HydPy| :class:`~hydpy.core.hydpytools.HydPy`
        """
        if self.master is not None:
            self.master._medium2long.update(self._medium2long)
            self.master.update_masters()

    def get_commands(self, source=None):
        """Return a string containing multiple `reStructuredText`
        replacements with the substitutions currently defined.

        Some examples based on the subpackage |optiontools|:

        >>> from hydpy.core.autodoctools import Substituter
        >>> substituter = Substituter()
        >>> from hydpy.core import optiontools
        >>> substituter.add_module(optiontools)

        When calling |Substituter.get_commands| with the `source`
        argument, the complete `short2long` and `medium2long` mappings
        are translated into replacement commands (only a few of them
        are shown):

        >>> print(substituter.get_commands())
        .. |Options.checkseries| replace:: \
:const:`~hydpy.core.optiontools.Options.checkseries`
        .. |Options.dirverbose| replace:: \
:const:`~hydpy.core.optiontools.Options.dirverbose`
        ...
        .. |optiontools.Options.warntrim| replace:: \
:const:`~hydpy.core.optiontools.Options.warntrim`
        .. |optiontools.Options| replace:: \
:class:`~hydpy.core.optiontools.Options`

        Through passing a string (usually the source code of a file
        to be documented), only the replacement commands relevant for
        this string are translated:

        >>> from hydpy.core import objecttools
        >>> import inspect
        >>> source = inspect.getsource(objecttools)
        >>> print(substituter.get_commands(source))
        .. |Options.reprdigits| replace:: \
:const:`~hydpy.core.optiontools.Options.reprdigits`
        """
        commands = []
        for key, value in self:
            if (source is None) or (key in source):
                commands.append('.. %s replace:: %s' % (key, value))
        return '\n'.join(commands)

    def find(self, text):
        """Print all substitutions that include the given text string."""
        for key, value in self:
            if (text in key) or (text in value):
                print(key, value)

    def __iter__(self):
        for item in sorted(self._short2long.items()):
            yield item
        for item in sorted(self._medium2long.items()):
            yield item


@make_autodoc_optional
def prepare_mainsubstituter():
    """Prepare and return a |Substituter| object for the main `__init__`
    file of :ref:`HydPy`."""
    substituter = Substituter()
    for module in (builtins, numpy, datetime, unittest, doctest, inspect, io):
        substituter.add_module(module)
    for subpackage in (auxs, core, cythons):
        for dummy, name, dummy in pkgutil.walk_packages(subpackage.__path__):
            full_name = subpackage.__name__ + '.' + name
            substituter.add_module(importlib.import_module(full_name))
    substituter.add_modules(models)
    for cymodule in (annutils, smoothutils, pointerutils):
        substituter.add_module(cymodule, cython=True)
    substituter._short2long['|pub|'] = ':mod:`~hydpy.pub`'
    substituter._short2long['|config|'] = ':mod:`~hydpy.config`'
    return substituter


def _number_of_line(member_tuple):
    """Try to return the number of the first line of the definition of a
    member of a module."""
    member = member_tuple[1]
    try:
        return member.__code__.co_firstlineno
    except AttributeError:
        pass
    try:
        return inspect.findsource(member)[1]
    except BaseException:
        pass
    for value in vars(member).values():
        try:
            return value.__code__.co_firstlineno
        except AttributeError:
            pass
    return 0


@make_autodoc_optional
def autodoc_module():
    """Add a short summary of all implemented members to a modules docstring.

    Just write `autodoctools.autodoc_module()` at the very bottom of the
    module.

    Note that function |autodoc_module| is not thought to be used for
    modules defining models.  For base models, see function
    |autodoc_basemodel| instead.
    """
    module = inspect.getmodule(inspect.currentframe().f_back.f_back)
    doc = module.__doc__
    if doc is None:
        doc = ''
    lines = ['\n\nModule :mod:`~%s` implements the following members:\n'
             % module.__name__]
    members = []
    for (name, member) in inspect.getmembers(module):
        if ((not name.startswith('_')) and
                (inspect.getmodule(member) is module)):
            members.append((name, member))
    members = sorted(members, key=_number_of_line)
    for (name, member) in members:
        if inspect.isfunction(member):
            type_ = 'func'
        elif inspect.isclass(member):
            type_ = 'class'
        else:
            type_ = 'obj'
        lines.append('      * :%s:`~%s` %s'
                     % (type_, name, description(member)))
    module.__doc__ = doc + '\n\n' + '\n'.join(lines) + '\n\n' + 80*'_'


autodoc_module()